Africanized honeybee population (Apis mellifera L.) in Nicaragua: Forewing length and mitotype lineages.

Various subspecies of Apis mellifera L. were introduced to Central America since colonization 500 years ago. Hybridization increased with the entrance of the Africanized bee in Nicaragua in 1984. Rustic beekeeping activities and numerous feral swarms define the genetic pattern, reflected in phenotypic heterogeneity and remarkable differences in the behaviour of the bee colonies, especially the nest defence. Due to these facts, the question emerge about the degree of Africanization of honeybee colonies in Nicaragua. In this study, we identified Africanized honeybee colonies based on the single character "mean forewing length" and we corroborated our results by determining mitotypes using mtDNA analysis. Morphometric and genetic approaches were realized in three different geographical zones of Nicaragua and related to beehive characteristics and management. Worker bee samples were taken from the inside of 146 hives from 26 apiaries. Abdominal colour as phenotypic character was the first examination, followed by measurement of 1460 right forewings to determine corresponding probability of Africanization. More than 60% of the beehives showed phenotypic heterogeneity and mean forewing length of 8.74 mm (SD 0.16 mm) indicated a high degree of Africanization. Those results provided a selection of 96 worker bees to perform PCR of two worker bees per hive. For mitochondrial DNA analysis 14 samples from sentinel apiaries were added. Three from 61 beehives presented bees with different mtDNA. Throughout, three mitotypes of the African (A) lineage were detected; one mitotype is still unidentified. Mitotype A1 A. mellifera iberiensis was represented by 88 bees and mitotype A4 A. mellifera scutellata by 21 bees. Phylogenetic analysis confirmed the PCR findings. No associations were found between mitotypes, forewing length, beehive characteristics and management. A high degree of Africanization in A. mellifera colonies represented by two predominating mitotypes from the A lineage, prevail in Neotropical Nicaragua, with mitotype A4 predominating at higher altitudes.

Introduction

The native range of the Western honeybee (Apis mellifera) is broad, comprising Africa, most of Europe and the Middle East. However, the evolutionary origin of the subspecies is still debated, whether it was in Middle East [1-4] or in Africa [5, 6] or even Middle East/North-eastern Africa [7].

Studies based on morphometric [3] and genetic analyses [8, 9], identified five evolutionary lineages of A. mellifera (hereinafter referred to as lineage A, M, C. O and Y). The identification of existing subspecies of the Western honeybee over the past 30 years has led to diverse results, due to the use of different identification approaches. With the Morphometric Bee Data Bank in Oberursel, Germany, Ruttner has created a solid framework for the classification of the subspecies of A. mellifera [3]. Up to 40 distinct morphometric characters, manually measured, were correlated to geographic variability for honeybee subspecies. Technical progress made computer assisted measurement and analysis possible [10-12], and ongoing morphometric research about A. mellifera species focused on geometric morphometric techniques using landmark-based features for subspecies identification [13, 14]. The combination of feature selection techniques (relevant morphometric characters for classification) and classifiers (machine learning algorithm) is used to improve the honeybee subspecies identification [15].

Genetic research offered new possibilities to infer phylogenetic relationships between honeybee species; in particular, mitochondrial DNA (mtDNA) analysis is effective due to maternal inheritance in bee colonies [1]. The intergenic COI-COII region of the mtDNA of A. mellifera is preferably investigated due to the variable number of a 192–196 bp sequence (Q) and the complete or partial deletion of the 67 bp sequence (Po). Length variability combined with a restriction site polymorphism enables a quick and easy test to characterize mtDNA mitotypes. Beside mtDNA analysis, nuclear markers are used in population genetic studies. The analysis of the microsatellites (simple sequence repeats, SSR) is able to identify the variation in the number and type of repetitions of the loci and thus can determine the genetic diversity in honeybee populations [8]. Single Nucleotide Polymorphism (SNP) are molecular genetic markers used for genetic diversity research. SNP analysis in bees, identify the sites where they differ in a single nucleotide in nuclear DNA sequence, leading to different genotypes [4]. In countries with lower budget for investigation, an alternative is the sequencing of the genetic code of the mtDNA intergenic COI-COII region [16]. In addition to the studies focused on one approach, investigations with morphometric and genetic approaches were conducted in order to identify subspecies of honeybees, to compare populations at different times and in different regions and to determine the process of Africanization [17-19].

Currently, the number and geographical distribution of honeybee subspecies is not clear due to population dynamics in transition zones [20] and the human impact on genetic diversity due to commercial activities [21], nonetheless the present taxonomic classification of A. mellifera mainly describes ~ 30 subspecies [5, 21].

In America, the evolutionary process did not take place as it did in Europe and Africa, because A. mellifera is not originated there. Honeybees came to America with European colonists and have spread across the continent for the past 500 years. The European subspecies brought to America were likely A. mellifera mellifera (France, Germany) and A. mellifera iberica (actually A. m iberiensis) and later the Italian bees, A. mellifera ligustica [22].

If the introduction of the Western honeybee to the American continent by the European colonialists was a successful biological invasion, the later introduction of the African subspecies into Brazil became an even more decisive event. In 1956, tropical African queens of the subspecies A. mellifera scutellata were introduced to exchange genetics with resident European stocks of honeybees to improve honey production [23]. By crossing European honeybee species with African queens and drones, the Africanized bee emerged and then spread into South, Central and North America [13, 24].

Quality methods measuring a significant number of morphometric characters have discriminate power to identify subspecies, but they are time consuming. Meixner et al. outlined the documented methods for discrimination of morphometric subspecies and lineages of the honeybee by measuring size characters of different body parts, colouration patterns, characteristics of pilosity and a detailed wing shape analysis [19]. In order to differentiate between Africanized and European honeybees as a study purpose, Rinderer et al. developed a simplified method by measuring the forewing length. This single character compared to multivariate analysis showed to be the best criterion to identify Africanization (at p<0.1 level no misidentification) with the overall average of the forewing length of an Africanized bee colony = 8.87 mm and = 9.20 for European bees [25]. Even assuming that this technique provides only preliminary results [19] and that it is unsuitable for identifying subspecies, it does indicate the probability that a colony is Africanized [26].

Genetic investigations are more accurate to discriminate the African and European matriline identities using PCR-RFLP and sequencing analysis. Widely used is the study of the intergenic region of mitochondrial DNA (mtDNA) COI and COII, because of its high degree of genetic variation between the different lineages. The region is composed of a P and/or Q sequence; depending on presence or absence (P-sequence), as well numbers and variation (Q-sequence) the mitotypes within each lineage can be determined [27, 28]. Studies with different restriction enzymes for digestion were conducted during the past 3 decades [18, 19, 29], but the DraI method showed a significant amount of polymorphism detection and it is highly suitable to discriminate Africanized colonies [27, 28].

In Nicaragua, the Africanized honeybee was first reported in 1984 [30]. Due to the high defensiveness of the Africanized bees, in the following years beekeeping activity was reduced over 50% (C. Düttmann personal communication to beekeeper committee), like in many other Latin-American countries where the Africanization of A. mellifera had occurred [18]. Since Nicaraguan beekeepers do not practice selective breeding and the majority catch feral swarms to increase the number of their hives, a process of natural selection began for the benefit of the Africanized bee, which replaced the European colonies nationwide. Despite the high level of defensiveness and swarming behaviour, the beekeepers got to know the advantages of working with the Africanized bee, for example resistance and tolerance to diseases, high reproductive numbers and higher honey yields under appropriate management conditions [30, 31].

Although the presence of the Africanized honeybee in Nicaragua is well known, there is still no published information about their morphology and genetic variability. It is important to study the genetic structure of the honeybee populations in Nicaragua, because this information could be used in future genetic improvement programs including the positive traits and attenuating the negative ones of the Africanized honeybee.

The present study wants to give a first insight into the Nicaraguan honeybee population that due to its behaviour characteristics probably has a high Africanized level. The applied methods are measuring of the forewing length, and the determination of mtDNA by PCR and sequencing.

We hypothesize that honeybee population from Nicaragua are mainly composed of Africanized bees, due to rustic apiculture management, where beekeepers capture feral swarms in order to introduce these colonies into their apiaries and furthermore, the practice to exchange frequently combs with brood from strong hives to reinforce weak hives [26, 32, 33].

Methods

Study sites

Nicaraguan beekeeping activities are dominated by a tropical savanna climate, with a distinction between three regions characterized by different temperature and precipitation averages and variation in altitudes of the apiary sites. The Pacific Lowlands with a mean temperature of 26.8°C reaching 45°C as maximum temperature (tropical savanna climate) and with apiaries at low elevations sites from 20 to 100 meters above sea level (m), exceptionally apiaries on the volcanic chain at 290 m. The Central Zone with lower to similar temperatures (mean 24.9°C), but tropical monsoon climate during the summer time, and higher elevations than the Pacific Lowlands with apiaries sites from 165 to 375 m. The Northern Highlands with the lowest temperatures (mean 20.4 to 22.8°C with minimum averages of 14°C) also tropical monsoon climate and apiaries locations from 280 up to 1009 m [34].

Sample collection

In order to obtain the specimens for the subspecies identification of A. mellifera in three geographically different zones of Nicaragua, samples of worker bees were taken from 146 beehives (26 apiaries): 67 beehives (7 apiaries) in the Pacific Lowlands, 38 beehives (6 apiaries) in the Central Zone and 41 beehives (13 apiaries) in the Northern Highlands. There is more beekeeping activity in the Pacific Lowlands, but the northern region has an extensive variety of microclimates and flora providing different conditions for beekeeping. Samples were collected between 2013 and 2016.

The sample of ~ 250 worker bees per colony was taken from three different frames of the breeding chamber, stored in labeled plastic containers (500 ml) containing 150 ml 96% ethanol and taken to the laboratory of the Centro Veterinario de Diagnóstico e Investigación, (CEVEDI), Escuela de Ciencias Agrarias y Veterinarias (ECAV), Universidad Nacional Autónoma de Nicaragua-León (UNAN-León).

Data collection about defensive behaviour and beehive management

During the sampling, additional data were registered, such as the defensive behaviour of the colony measured by massive attack during the sampling and persecution after leaving the apiary indicated in meters (modified method, defensive response [35]. The ordinal classification of the defensive behaviour: 1. gentle: no massive attack and low persecution (<30 m) by few bees; 2. moderate: no massive attack, but strong persecution (>30 m) or massive attack, but low persecution (<30 m); 3. very defensive: immediate and high massive attack, followed by strong persecution (>30 m up to more than 100 m). Management information was recorded from the beekeepers about the applied methods of multiplication of the beehives, queen bee replacement and the origin of the queen in the beehive. To collect data, the investigators developed a questionnaire that was pilot–tested before application. The flow chart of investigation procedure is shown in Fig 1.

Fig 1

Flow chart of investigation procedure.

Part 1: Phenotypic and morphometric identification

Although it is not used as a main discrimination factor, the abdominal colour variability still is a visual sign for intracolonial diversity, distinguished by the pigmentation variation of the abdominal tergites and the scutellum. The simple observation of the bimodal variation between “yellow” and “dark” honeybees in the same beehive indicates phenotypical diversity. As the first feature in this study, phenotypic diversity of the bee colonies in each of the 146 hives was examined using photographic identification of the collected bees to document uniformity or diversity of the abdominal colour pattern. [3, 19].

As an effective screening tool, the measurement of the length of the right forewing was used as a single character and simple morphometric identification that best discriminate Africanized and European honeybees, as described by Rinderer et al. (1986). From each sample, ten specimens were selected; in the case of phenotypic heterogeneity, half of the bees were dark and the other half yellow. The right forewing of each selected individual was removed as close to the thorax as possible and mounted on a slide, stereoscopically photographed, digitalized and measured with the program Image Tool 3, (S1 Fig), [36].

Part 2: Molecular mtDNA identification

For molecular diagnosis, honeybees were selected from all inspected apiaries in study. In order to obtain relatively high genetic diversity from the Nicaraguan honeybee population, the selection criteria based on heterogeneity of abdominal colour within the colony and mean length of forewing. Due to the management of beekeepers in Nicaragua, referring to frequently exchanging combs with brood to strengthen weak hives, exists the probability of finding diversity in mitochondrial DNA in one beehive [30]. Therefore, we decided to sample two worker bees per beehive. Phenotypically there is a high diversity within the colonies and in the recollected samples, which might not be only because of polyandry. As Evans et al. described, that in case selected worker bees could be from different colonies, it is worth sampling more than one individual to avoid mistakes in assigning colony heritage [22]. The decision to select more than one worker bee is a modification of the usual methodology in processing mtDNA with the DraI-Test.

Additionally, only for genetic analyses (but not for forewing length and additional data), honeybees from five sentinel apiaries of the annual monitoring program from the Instituto de Protección y Sanidad Agropecuaria (IPSA) were included, because of the particular location of the apiaries: northern border to Honduras (1), southern border to Costa Rica (1), Island in the Nicaraguan Lake (1) and the South Caribbean Coast (2). Inclusively, the sample size for PCR selection comprised 110 worker bees from 61 beehives from 31 apiaries; and per zone Pacific Lowlands (33), Central Zone (26) and Northern Highlands (51).

DNA extraction

The 110 bees preserved in 96% ethanol and stored at −20°C until molecular analysis, were washed three times with PBS and once with distilled water. Total DNA isolation was performed from the two hind legs of each individual following the method of Garnery et al. (1993) and according to the protocol described in the Qiagen Blood and Tissue DNA kits (Qiagen, Germany)) with an adjusted elution volume of 50 μl of each sample instead of 200 μl.

Mitochondrial DNA analysis

For PCR amplification of the intergenic COI-COII region E2 forward primer 5’-GGCAGAATAAGTGCATTG-3’ and H2 reverse primer 5’-CAATATCATTGATGACC-3’ were used. The reaction was carried out in a total volume of 50 μl, and contained 25 μl of Master Mix 2X (Promega, USA), 2 μl each primer (10,000 nM), 16 μL of Nuclease-free water and 5 μl of DNA of the samples. The amplification protocol consisted of an initial temperature of 95°C for 5 minutes, followed by 35 cycles of: 95°C for 30 seconds, 47°C for 30 seconds and 72°C for 1 minute. Subsequently a final elongation at 72°C for 15 minutes and hold at 4°C. The PCR reaction was performed with the Applied Biosystem 2720 Thermocycler. The amplification products were separated by placing 10 μl of the mixture on 1.5% agarose gel with ethidium bromide and visualized under UV light.

The PCR products were digested with DraI restriction enzyme to identify mitotypes of A. mellifera at the recognition and cleavage sequence TTTAAA. According to the protocol from Thermo Fisher Scientific (MA, USA), adding 18 μl of nuclease-free water, 2 μl of 10X buffer, 2 μl of the enzyme, 10 μl of the PCR product. It was mixed gently and incubated at 37°C for 16 hours. The digestion products were observed on a 1% agarose gel with ethidium bromide. The determination of the mitotypes was carried out according to the number and size of the DNA segments as described by Garnery et al. (1993).

Part 3: Sequencing

To confirm the findings in PCR diagnosis, ten specimens were selected for sequencing from the three sample zones: three specimens from the Pacific Lowlands, three from the Central Zone and four from the Northern Highlands. In addition to their geographical origin, the ten specimens were selected based on the results of the determination of the mitotype; likewise, it was considered that in the Central Zone and Northern Highlands the result of the PCR indicated different mitotypes in the same hive.

Reverse-sense sequencing (H2 primer) was performed for the ten selected samples. DNA purifications were carried out with the commercial kit according to the manufacturer UltraClean® 15 DNA Purification (MO BIO, USA). The sequencing was carried out at the Centro de Investigación en Biología Celular y Molecular (CIBCM), Universidad de Costa Rica (UCR).

Phylogenetic analysis

The sequences of each sample were aligned using ClustalW 1.6 to study evolutionary history. Evolutionary divergences were calculated using the Maximum Likelihood (ML) method and reported as units of the number of base substitutions per site. Bootstrap test (1000 repetitions) was performed to determine statistical significance and expressed as the percentage of occurrence of each taxon in the branches of the phylogenetic tree.

Evolutionary analysis was conducted in MEGA 7 software [37]. The sequences were sent to the National Center for Biotechnology Information (NCBI) using Banklt (NCBI, 2021), obtaining the identification numbers: MW695396, MW695397, MW695398, MW695399, MW695400, MW695401, MW695402, MW695403, MW695404, and MW695405.

Statistical analysis

Qualitative variables were expressed in relative frequency. ANOVA test was realized for the mean forewing length related to zones, departments, municipalities of sampling and characteristics of beehives and their management. The degree of Africanization of the beehives was determined according to the table “Fore-wing length means ( European (EA)” [25]. Pearson correlation test was applied to relate forewing length means with altitude of sampling sites.

Student’s t test was applied to compare quantitative variables with normal distribution such as the length of the forewing and the altitude of the locations of the apiaries, both related to the two identified mitotypes.

Chi Square tests (χ2) were performed to show association between the geographic distribution of the mitotypes and characteristics of the beehives and their management.

A Fisher’s exact test was applied to compare the frequency of mitotypes between ten sampled departments of the country, and considered statistically significant at α = 0.05. Data management and analyses were performed using R v.3.4.1 statistical software (Foundation for Statistical Computing, Vienna, Austria).

Ethical approval

The Research Ethics Committee (ECAV, UNAN-León) has previously approved the study. The samples were sourced from 26 private apiaries with the permission of the apiary owners; on the other hand, samples were provided from IPSA from epidemiological monitoring. The animals collected are not considered endangered or protected by any institution.

Results

The present study reveals a first insight to the Africanization level of honeybee colonies in Nicaragua by measuring the forewing length and identifying the different mitotypes in three zones with beekeeping activities, sustained by findings of the defensive behaviour. The reproduction management of the apiaries is a main activity that contributes to colony diversity, reflected by phenotypic diversity identified by abdominal colour pattern.

Colony characteristics and beehive management

The defensive behaviour of the studied 146 colonies was related to Africanized characters in A. mellifera. Three quarters of the colonies had a defensive tendency; 20.5% were highly defensive with massive attacks during handling the hives and persecutions from 30 up to 400 meters leaving the apiary; and 54.1% of the beehives were moderately defensive, requiring special management to avoid problems. Only 25.3% were gentle bees, without any problem in management.

The obtained data revealed four methods used for multiplication of the 146 beehives. The non-technical method was capturing feral swarms (22.6% of the hives) and 93.9% of them remained without queen replacement. Likewise, beekeepers bought hives from other beekeepers without replacement of the queen bee (17.8% of the hives). The most applied multiplication method was the division of the hives (46.6%); nearly all of the new introduced queens were self-reared by the beekeepers (98.5%).

In 57.5% of the studied hives, the beekeepers replaced the queens and 95.2% of those beekeepers used self-reared queens. Only 4.8% of the queens were purchased from one cooperative with special breeding techniques. There is no queen-rearing program in Nicaragua that could supply selected queens to beekeepers.

The hive characteristics and their management were significantly different in between the sampling zones (Table 1).

Table 1

Association of Apis mellifera hive characteristics and management according to sampled zones.

Hive characteristics and management		Sampled zones			DF	χ2	p
		Pacific	Central	North
Abdominal colour pattern	uniform	4	8	26	2	12.39	< 0.01
	different	25	11	22
Defensiveness	gentle	12	0	5	4	29.06	< 0.01
	moderate	13	17	21
	very defensive	4	2	22
Origin of the beehive	captured	14	0	3	6	89.32	< 0.01
	bought	3	9	4
	divided	12	1	41
	reared	0	9	0
Queen replacement	yes	12	10	38	2	12.00	< 0.01
	no	17	9	10
Origin of the queen	feral	17	0	6	4	35.80	< 0.01
	own	12	17	42
	reared	0	2	0

DF: Degree of freedom, χ2: Statistic value, p: value

Phenotypic diversity was detected in 62.3% of the hives according to the pattern of abdominal colour with variations from dark brown to light yellow in the same colony (Fig 2); there was a higher degree of heterogeneity among the hives in the Pacific Lowlands (Table 1).

Fig 2

Worker bees from one beehive (Apis mellifera) in Boaco, Central Nicaragua.

Seven worker bees from the same sample of one beehive showing high polymorphism according to abdominal colour pattern.

Forewing length

The overall average of the right forewing length ( = 8.74 mm; SD = 0.16 mm), measured in 146 sampled beehives, indicated a high degree of Africanization. This study revealed that in Nicaragua more than 79% of the beehives had a degree of 99% of Africanization. The level of 90% of Africanization was detected in 87% of the beehives; there were only 8/146 managed colonies with a European character more than 50%, according to the forewing length (Fig 3). The mean length of the right forewing per beehive was not different among the sample sites (ANOVA, zone: F2,145 = 1.77, p > 0.05; department: F6,145 = 1.32, p > 0.05; municipality: F13,145 = 1.25, p > 0.05) (Table 2), also there was no correlation between forewing length and altitudes of the localization of the apiaries (R = -0.02, p > 0.05). In addition, there was no significant difference between the characteristics of the beehives and management activities according to the forewing length (ANOVA, origin of the colony: F3,146 = 0.13, p > 0.05; abdominal colour pattern: F1,146 = 0.55, p > 0.05; defensiveness of the beehive: F2,146 = 1.52, p > 0.05; queen replacement: F 1,146 = 0.48, p > 0.05; origin of the queen: F2,146 = 0.66, p > 0.05) (Fig 4).

Fig 3

Degree of Africanization in beehives (Apis mellifera) in Nicaragua.

Yellow dotted lines 90% of Africanization; Red dotted lines 99% of Africanization at beehive level.

Table 2

Beehive mean length of the right forewing from Apis mellifera at different location level.

Location level		sampled in		Mean Wing Length	SD	95% Confidence Interval		DF	F	p
		# Apiaries	# Beehives			lower bound	upper bound
Zone	Pacific	7	67	8.74	0.16	8.70	8.78	2	1.768	0.17
	Central	6	38	8.71	0.16	8.66	8.77
	North	13	41	8.78	0.17	8.73	8.83
Department	León	4	38	8.73	0.15	8.68	8.78	6	1.317	0.25
	Chinandega	3	29	8.75	0.16	8.69	8.81
	Boaco	6	38	8.71	0.16	8.66	8.77
	Matagalpa	5	15	8.78	0.21	8.66	8.90
	Jinotega	5	15	8.82	0.14	8.74	8.90
	Estelí	1	3	8.62	0.23	8.05	9.18
	Madriz	2	8	8.77	0.07	8.71	8.83
Municipality	León	3	28	8.72	0.15	8.66	8.78	13	1.250	0.25
	Nagarote	1	10	8.75	0.17	8.62	8.87
	Villanueva	3	29	8.75	0.16	8.69	8.81
	Boaco	1	5	8.71	0.13	8.54	8.87
	Teustepe	5	33	8.71	0.17	8.65	8.77
	Ciudad Dario	2	6	8.66	0.05	8.61	8.71
	Tuma—La Dalia	1	3	8.80	0.36	7.91	9.70
	Matiguás	1	3	8.81	0.23	8.23	9.40
	Muy Muy	1	3	8.94	0.21	8.42	9.47
	Jinotega	1	3	8.71	0.22	8.17	9.26
	El Cuá	3	9	8.85	0.13	8.75	8.95
	Wiwilí	1	3	8.84	0.06	8.69	8.99
	Estelí	1	3	8.62	0.23	8.05	9.18
	Somoto	2	8	8.77	0.07	8.71	8.83

SD: Standard deviation, DF: Degree of freedom, F: Statistic value, p: value significant at <0.05

Fig 4

Beehive characteristics according to sampled zone and related to forewing length.

(a) Forewing length by sampled zone; (b) Forewing length by origin of the colony and sampled zone; (c) Forewing length by abdominal colour pattern and sampled zone; (d) Forewing length by defensiveness of the hive and sampled zone; (e) Forewing length by queen bee replacement and sampled zone; (f) Forewing length by origin of the queen and sampled zone.

Mitotype identification

The results obtained by DraI test indicated that 100% of the studied A. mellifera population (110 worker bees) had an African mitotype, represented in three classes; 80% (88/110) of the samples were identified as mitotype A1, while 19.1% (21/110) were identified as mitotype A4. One mitotype is still unidentified and will be submitted to further identification. The distribution of the mitotypes by the sampled zones was significantly different (sampled zones: χ24,110 = 10.48, p<0.05). According to χ2, the mitotype A1 was found with a frequency of 80.8% (21/26) and 97% (32/33) in the Central Zone and Pacific Lowlands respectively, while in the Northern Highlands the frequency of this mitotype was 68.3% (35/51). Concerning the altitude of the apiaries, the distribution of the mitotypes was significantly different (p < 0.001), finding mitotype A4 mostly in higher locations ( = 674.78 m, SD = 281.14 m) and the more frequent mitotype A1 in lower zones ( = 316.68 m, SD = 275.76 m) (Fig 5).

Fig 5

Geographic distribution of Apis mellifera mitotypes by department, Nicaragua.

The association of the mitotypes with the characteristics of the beehives and their management was calculated based on the genetic results of the 96 bees from the main database; 14 mitotypes were excluded, because the information wasn´t complete (honeybee samples with courtesy of epidemiological surveillance by IPSA) (Fig 1). The analyses revealed, that there is no significant difference between the mean forewing length of the identified mitotypes (p > 0.05), mitotype A1 with a mean forewing length = 8.75 mm (SD = 0.18 mm) and mitotype A4 with = 8.79 mm (SD = 0.21 mm). No significant difference was observed linking mitotypes to the defensiveness of the hives and the phenotypic diversity (defensiveness: χ24,96 = 3.39 p≥0.05; abdominal colour pattern: χ22,96 = 2.68 p≥0.05). Neither did the comparison of the mitotypes with the management of the hives (origin hive: χ26,96 = 8.60 p≥0.05; replacement queen: χ22,96 = 0.80 p≥0.05; origin queen: χ24,96 = 3.11 p≥0.05).

Worker bees with different mtDNA in the same beehive were detected in 3/61 beehives, which belonged to one apiary in the Central Zone and two apiaries in the Northern Highlands.

Phylogenetic analysis

The phylogenetic tree constructed by the ML method based on inference analysis of the COI-COII intergenic region, classified three bees of the present study as A. mellifera scutellata mitotypes A4 (MW695396, MW695397, MW695398) in the same group other sequences were described from Africa (MG592308, FJ477987, KF977009, AF503562), Brazil (EF033650) and Argentina (JQ582439). The sequences of six samples were classified as A. mellifera iberiensis mitotypes A1 (MW695399, MW695400, MW695401, MW695402, MW695404, MW695405) clustered with samples from Africa (KX463745, MG592298), Europe (KX463736, FJ477985) and America (EF033649 JQ582438, FJ743639). The sequences corresponding to the mitotypes M4 (EF033656), C22 (FJ037785) and C1 (EF033655) were grouped distantly from the Nicaraguan mitotypes A4 and A1 (Fig 6).

Fig 6

Phylogenetic analysis of COI-COII gene nucleotide sequences of Apis mellifera.

The evolutionary history was inferred by using the Maximum Likelihood method and Tamura-Nei model [38]. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The proportion of sites, where at least one unambiguous base is present in at least one sequence for each descendent clade, is shown next to each internal node in the tree. This analysis involved 29 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. There were a total of 656 positions in the final dataset.

The evolutionary divergence between the three samples from Nicaragua classified as mitotype A4 ranged from 0.000 to 0.003 (similarity more than 99.7%), same result with the sequences from other countries with similarity more than 99.7, exceptionally one Nicaraguan mitotype A4 (MW695396) that was 100% similar to sequences described in South Africa, Brazil and Argentina (FJ477987, EF033650, JQ582439). The similarity between the six samples from Nicaragua classified as mitotype A1 was at least 99.8%; five of these samples were 100% similar to the A1 mitotype described in Argentina (JQ582438), and four were 100% similar to the sequence from Benin (MG592298). Relating A1 mitotypes from south-western Europe to the Nicaraguan A1 sequences, the divergence was higher (between 0.002 and 0.004), even higher the comparison with the sequence from North Africa (between 0.004 and 0.006). The sample MW695396 (Mitotype A4) and the sample MW695399 (Mitotype A1) were collected from the same hive and presented a divergence of 0.119 (Table 3).

Table 3

Estimation of evolutionary divergence between selected Apis mellifera mitotypes.

		1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19
1	A4-MW695397-Nicaragua
2	A4-MW695396-Nicaragua	0.003
3	A4-MW695398-Nicaragua	0.000	0.003
4	A4-MG592308-Benin	0.002	0.002	0.002
5	A4-FJ477987-South Africa	0.002	0.000	0.002	0.000
6	A4-EF033650-Brazil	0.003	0.000	0.003	0.002	0.000
7	A4-JQ582439-Argentina	0.003	0.000	0.003	0.002	0.000	0.000
8	A1-MW695405-Nicaragua	0.118	0.108	0.118	0.112	0.108	0.108	0.108
9	A1-MW695404-Nicaragua	0.117	0.107	0.117	0.111	0.107	0.107	0.107	0.000
10	A1-MW695402-Nicaragua	0.117	0.107	0.117	0.111	0.107	0.107	0.107	0.000	0.000
11	A1-MW695401-Nicaragua	0.116	0.110	0.116	0.114	0.111	0.110	0.110	0.002	0.002	0.002
12	A1-MW695400-Nicaragua	0.117	0.107	0.117	0.111	0.107	0.107	0.107	0.000	0.000	0.000	0.002
13	A1-MW695399-Nicaragua	0.119	0.110	0.119	0.114	0.110	0.110	0.110	0.000	0.000	0.000	0.002	0.000
14	A1-EF033649-Brazil	0.193	0.190	0.210	0.206	0.205	0.206	0.208	0.002	0.002	0.002	0.004	0.002	0.002
15	A1-JQ582438-Argentina	0.195	0.188	0.195	0.191	0.188	0.188	0.188	0.000	0.000	0.000	0.002	0.000	0.000	0.003
16	A1n-KX463745-Algeria	0.199	0.196	0.216	0.212	0.211	0.212	0.214	0.004	0.004	0.004	0.006	0.004	0.004	0.003	0.005
17	A1e-MG592298-Benin	0.197	0.190	0.213	0.209	0.207	0.207	0.208	0.000	0.000	0.000	0.004	0.000	0.002	0.005	0.003	0.005
18	A1b-FJ477985-Spain	0.168	0.164	0.168	0.163	0.162	0.163	0.164	0.002	0.002	0.002	0.004	0.002	0.002	0.000	0.002	0.002	0.003
19	A1-KX463736-Portugal	0.197	0.194	0.214	0.210	0.209	0.210	0.212	0.002	0.002	0.002	0.004	0.002	0.002	0.002	0.003	0.002	0.003	0.000

The number of base differences per site from between sequences are shown. This analysis involved 19 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All ambiguous positions were removed for each sequence pair (pairwise deletion option). There were a total of 656 positions in the final dataset.

Discussion

The process of the Africanization of the honeybee A. mellifera in the continent of the Americas since the 60´s of the last century began in Nicaragua more than 3 decades ago. The abundant flora of the country and a mostly rustic beekeeping with little intervention from beekeepers, helped honeybees with African characteristics to reproduce easily, especially, when colonies of feral bee populations formed part of the stock managed by beekeepers [18, 26]. Table 1 shows that there are significant differences between the sample zones in terms of defensiveness, the origin of the hive, replacement of the queen and the origin of the same. Curiously, the results reflect that in the Pacific Lowlands bees are docile, although the capture of feral swarms is widely practiced. In contrast, in the Northern Highlands offensive behaviour predominates, although the division of hives with queen replacement is common. It appears that defensive colonies are being reproduced, a finding that may in part be explained by the interest of beekeepers in avoiding theft of their hives. Queen rearing is merely practiced in the Central Zone of the country; also, there are few beekeepers who buy queens to replace them in new hives. As mentioned before, in the Pacific Lowlands, feral swarms are captured to multiply the hives, this is also reflected in the phenotypic heterogeneity. For this reason, it was decided to select two worker bees per sample for the identification of the mitotypes [27].

Due to these conditions and the results of the present study, in which A. mellifera from lineage A coincide with the morphometric character of the mean forewing length = 8.74 mm, a high degree of Africanization can be found in colonies of A. mellifera in Nicaragua, confirmed by the discrimination of Rinderer et al. (1986) who described = 8.87 mm as the overall average of the forewing length for Africanized beehives.

As Nicaragua has tropical climate, even in higher locations, it is not surprising to find Africanized honeybees across the country. There was no morphological variation nationwide, considering the mean forewing length as an indicator for Africanization [25], since there was no difference between the study areas or the height of the location of the apiaries (Table 2), neither association to characteristics of the hives and nor to their management (Fig 4). Amakpe et al. in Benin, Africa, found no correlation between morphometric characters (mainly leg and wing length) and altitude of the inspected apiaries, although they described a positive south-north gradient in bee size. Moreover, they indicated that there was no difference between the morphometric characters between the different haplotype, same results as in the present study in Nicaragua [39].

This study revealed that all the analysed bees belong to the A lineage, whereas 80% (88/110) of the samples were identified as mitotype A1, and 19.1% (21/110) as mitotype A4. In one sample, mitotype A could not be classified and will be submitted for further identification. There was a significant difference between mitotypes and altitudes corresponding to different zones, more honeybees with mitotype A4 were detected in higher altitudes ( = 674.78 m) of the Northern Highlands than in the Pacific Lowlands and Central Zone (Fig 5). Similarly, Kraus et al. described mitotype A4 more frequently in higher altitudes in Veracruz, Mexico. As a reference of an African population, they cited Franck et al. (2001) reporting 70% of mitotype A4 and 20% A1 at an altitude of 1300 m in Pretoria, South Africa with humid subtropical climate according to Köppen Geiger Climate Classification [40].

In Costa Rica with environmental conditions similar to those of Nicaragua, Lobo found a higher frequency of mitotype A1 (55.9%) than mitotype A4 (33.9%) [41]. In Colombia, 98.3% of the studied bee population had African mitotypes, represented in 17 classes (A), being A1 and A4 mitotype, the most frequent within this lineage [9]. In apiaries from Brazil and Uruguay, Collet et al. reported the predominance of the A1 and A4 mitotypes; however, they found a higher frequency of the A4 mitotype (68%) [42]. This difference can be attributed to the fact that generally A1 mitotype is more common in northern regions, both in Africa and Latin America, while the A4 mitotype has been observed more frequently in South Africa and South America [8, 42]. Further north from Nicaragua, in Mexico, Dominguez-Ayala et al, studied honeybee genetics in five main beekeeping regions with different climate. They reported colonies with mitotypes from lineages A, C and M; more than 50% of the colonies presented honeybees with mitotypes from the A lineage, mostly in subtropical and tropical regions. However, A1 was the most frequent mitotype in their study [17]. Likewise, Kraus et al. found in Veracruz, with tropical savanna climate in Mexico, honeybee mitotypes from the lineages A, C and M, but principally A1 and A4. In the lower part of the study site, samples were taken from hives established from locally captured swarms, indicating A mitotypes (67%) and European mitotypes (33%), meanwhile in the high altitudes (up to 2500 m) only feral colonies were sampled and the bee population was a 97% Africanized. They discussed these findings, that in the lower zones some beekeepers introduced constantly European queens, whereas in the higher altitudes Africanized bees replaced European bees, because of missing beekeeping activities [40].

The low diversity of mitotypes of the A lineage reported in this study could be explained by the action of the restriction enzyme DraI that is known to miss genetic variations in the A lineage [8, 9, 43] and still there is one unidentified mitotype from the Northern Highlands.

The sequencing of the intergenic region of the COI-COII confirmed the results of identified mitotypes by digestion with the restriction enzyme DraI. Both mitotypes A4 and A1 showed divergence between the samples from Nicaragua. Three bees in the present study were classified by sequencing as A. mellifera scutellata mitotype A4. The phylogenetic analysis indicated a similarity > 99.7% to sequences from Africa [39], Brazil [42] and Argentina [44], reinforcing the findings that this mitotype is found more frequently in subtropical and tropical areas [42].

Of the six mitotypes A1 from Nicaragua, five presented 100% similarity in between them, while one sample (MW695401) presented divergence (0.002). In this case, sequencing comparison is useful to identify genetic variants within the same mitotype, which cannot be evidenced by digestion pattern of the DraI test [39]. Remarkably, these five sequences were 100% similar to mitotype A1 from Argentina [45] and 4/5 to an A1 sequence from Benin, Africa [39]. Meanwhile the similarity to samples from south-western Europe and North Africa was lower (>99.6% and 99.4% respectively) [8, 46], not supporting that Nicaraguan A1 mitotypes may reflect Iberian honeybees introduced by the colonists [8]. In this discussion about the ‘pre-Africanization’ population, referring to mitotype A1 as A. mellifera iberiensis from the Iberian Peninsula, Kraus et al. described that in samples collected in Tamaulipas, Mexico, before the arrival of the Africanized honeybee in 1988, 13% of the honeybees presented mitotype A1 [40]. On the other hand, Chávez-Galarza et al. referred to “Ancestral Haplotypes” as mitotypes introduced from Africa to the Iberian Peninsula thousands of years ago, as part of the natural migration of A. mellifera [46]. Numerous studies about the subspecies from the Iberian Peninsula provide a complex dataset about phylogeographic patterns of A. mellifera [47], but still there are many questions about the origin of actual mitotypes from the A lineage detected in the Americas.

One characteristic for A1 mitotype belonging to A. mellifera iberiensis from south-western Europe could be the size of the forewing, commonly European subspecies have longer wings [48]. Conversely, in Nicaragua, colonies with mitotypes A4, generally characterized by smaller wings, here had even longer forewings than A1. These results cast doubt on the assumption that those bees are descendants from European honeybees, but only forewing length is not reliable to discriminate definitive subspecies, approaches like geometric morphometrics and if affordable combined to nuclear DNA markers, offer more information about the evolutionary processes [47].

Another important finding of this study disclosed, that sample MW695396 (mitotype A4) and sample MW695399 (mitotype A1) were collected from the same hive and presented a divergence of 0.119. In theory, all worker bees in a colony should have the same mtDNA, but beekeepers’ practices, like interchanging frames from brood chambers, may contribute to find more than one mitotype per hive. That means, sampling only one worker bee per hive may be insufficient to estimate the diversity of mitotypes [27].

Conclusion

In Nicaragua, regardless of the zone, there is a high degree of Africanization in A. mellifera colonies, represented by the mitotypes A1 and A4 and supported by the overall average of the right forewing length. The tropical climate prevailing throughout the country is a relevant factor for the habitat of Africanized bees, with mitotype A4 predominating at higher altitudes. At the same time, the colonies with the most defensive behaviour were reported in the Northern Highlands; defensiveness as a typical Africanized trait for honeybee colonies of the A4 lineage such as found in A. mellifera scutellata colonies. Beekeeping activities like queen replacements show a high impact in honeybee genetics, and this is probably why in Nicaragua, due to very rustic beekeeping, there has been a complete replacement of European bees by Africanized bees, at least no European lineages were detected by mtDNA analysis in the present study.

Nevertheless, Africanization in bee colonies presents advantages over European bees, related to resistance and/or tolerance to certain diseases mainly varroosis. These abilities are associated to a better hygienic behaviour and grooming, as well as a lower susceptibility to suffer invasion and reproduction of pathogens [9, 30, 31].

Measurement of the right forewing with image tool 3.

(TIFF)

Click here for additional data file.

Data information, Africanized honeybee population in Nicaragua, forewing length and mitotype lineages.

(XLSX)

Click here for additional data file.

10 Dec 2021

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Tiago M. Francoy, Ph.D.

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

3. We note that Figure 2 in your submission contain copyrighted images. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

a. You may seek permission from the original copyright holder of Figure 2 to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

b. If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

4. We note that Figure 5 in your submission contain [map/satellite] images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

a. You may seek permission from the original copyright holder of Figure 5 to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

b. If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

The following resources for replacing copyrighted map figures may be helpful:

USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/

The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/

Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.html

NASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/

Landsat: http://landsat.visibleearth.nasa.gov/

USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#

Natural Earth (public domain): http://www.naturalearthdata.com/

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: General comments:

The manuscript discusses the Africanization level of honeybee (Apis mellifera 1 L.) from Nicaragua using phenotypic, molecular, behavior, and morphometric analyses. I believe that this manuscript can be published in this Journal, but it is necessary to do a major reviewer. The authors need to organize better all methods that they used. They also need to work more in the discussion.

I suggest the authors add geometric morphometrics analysis to study shape and centroid size of the forewings because it can provide more accurate results than the traditional morphometrics that they used.

#Abstract

The authors need to work more on the abstract and include the question and the main objective of the study in the first paragraph.

Line 28: the authors used four methods to study honeybee populational variation. Please add the phenotypic and behavior in this line.

Line 28: Replace “approach” with approaches

# Introduction

Line 48: Replace “approach” with approaches

Line 51: replace “analysis” with analyses

Line 51-56: please extend this literature reviewer

Line 102: Choose “behaviour” (UK) or “behavior” (USA) in your manuscript.

Line 115-118: which previous studies did you use to elaborate your hypotheses? Please, include them.

#Methods

Line 133: are bees collected in the same year?

Lines 145-154: I suggest this paragraph be a subtitle about behavior

Lines 151-153: “Further information was recorded about the applied methods of multiplication of the beehives, queen bee replacement and the origin of the queen in the beehive”. What did you use to get this information?

Line 163: Are the “146 samples” the beehives? Please, specify it.

Line 166-172: Please provide an image of the forewing indicating the measurements. You also have two variables: colour and wing length. How many bees of the same colour did you use to measure the forewings? It is interesting to investigate if there are differences among bee forewings with different behaviors.

Line 172: Add the references of the software.

Line 176-182: why did you include samples in your work just for molecular analyses?

Lines 183-186: why did you include a new hypothesis here? Which previous studies did you use to elaborate this hypothesis?

Results

Line 263: you also studied the abdominal colour

Lines 274-283: the authors need to describe this procedure in the methods. How did you get this result?

Line 294: Write the species name in italic. Please, look at all manuscript.

Lines 300-306: present only your results and move this information to discussion.

Discussion

The authors need to discuss all result topics including the evolutionary history of Apis mellifera.

Lines 480-482: the authors agreed that geometric morphometrics can “offer more information about the evolutionary processes”. However, they used traditional morphometrics.

Reviewer #2: The process of the Africanization of Apis mellifera and its occupation in the Americas is a fascinating phenomenon that still needs to be studied. Information on the genetic and morphological diversity of the honeybee in Nicaragua contributes to a greater understanding of this process.

Here are some considerations:

- On line 55, it would be interesting to mention what these recent techniques are about. Unfortunately, in the text, it was not clear.

- Among a series of measures proposed by Ruttner in traditional morphometry to distinguish Apis mellifera lineages, do you think choosing only one linear measurement (affected by environmental factors) is enough to discriminate the groups? For example, according to Diniz-Filho et al., 2000 “Phylogenetically, the wing venation is more informative compared to the more environment–sensitive character categories of size, colour or pilosity”.

- In the text, the average length of the wing of European honeybees was not mentioned. Are they larger or smaller than the measurements of Africanized bees?

- Why was the entire sample of phenotypic variation not included in the molecular analyses? Despite the description of the criteria “...phenotypic diversity within the colony and mean length of the forewing.”, which factor determined this choice was unclear.

- On line 197, describe the modifications in the protocol.

- Neighbor-joining is a clustering algorithm that clusters haplotypes based on genetic distance. Maximum likelihood can apply a model of sequence evolution and is better for building a phylogeny using sequence data. Neighbor-joining isn’t a phylogenetic analysis

-Paragraphs 262-266 would fit better at the end of the text as a conclusion to the study.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

10 Jan 2022

León, 10th January 2022

Manuscript Number: PONE-D-21-31700

Dear Tiago M. Francoy, Ph.D.

Academic Editor

PLOS ONE

We are resubmitting our manuscript Africanized honeybee population (Apis mellifera L.) in Nicaragua: Forewing length and mitotype lineages. to PLOS ONE and all corrections are included. All the authors have reviewed the draft and approved its contents and conclusions.

Best regards,

The authors.

Dear Academic Editor, dear Reviewers.

Thank you for considering publishing our manuscript and for believing in the scientific merit of our contribution through this investigation. We have carefully revised your comments and suggestions and made adjustments accordingly to improve the manuscript. If you do not mind, we would first like to explain briefly the scientific situation of beekeeping in Nicaragua.

As a few academics in Veterinary Medicine from UNAN-León (National Autonomous University of León), we started studies on varroosis in 2004, following a personal request from some beekeepers. At the time, we were not beekeeping specialists and we had no previous references from published articles on beekeeping issues in Nicaragua. Therefore, we started to learn about these insects and their importance for us.

With the material and data from a cross-sectional study between 2013 and 2016, we had the first opportunity to contribute scientific information on beekeeping topics in Nicaragua. For the past 10 years, we worked together with specialists in Latin America to share our knowledge with the aim of improving our studies, above all to offer scientific solutions to practical problems in Nicaraguan beekeeping.

We have learned from “our bees” that they are very healthy and strong and that they have been able to defend themselves against diseases so far. (see article Düttmann C, Flores B, Sheleby-Elías J, Castillo G, Osejo H, Bermudez S, et al. Morphotype and haplotype identification of Varroa destructor (Acari: Varroidae), and its importance for apiculture in Nicaragua. Exp Appl Acarol. 2021;83: 527–544. doi:10.1007/s10493-021-00603-9).

In addition, we have learned from the Nicaraguan beekeepers that most of them practice rustic beekeeping, but as academic supervisors, we have good contact and acceptance with our observations on improving management and honey production. The process of Africanization in Nicaragua is of global concern. Rustic beekeeping mainly leaves it to nature to influence the population genetics of Apis mellifera. We therefore expect an almost “pure” Africanized bee population.

Your principal question why we didn´t applied geometric morphometrics is that in the first place we wanted to realize an estimation about the Africanization level in the honeybee population of Nicaragua, expecting a high degree. Therefor the method of Rinderer is the less time and money consuming one with valid results. On the other hand, we wanted to confirm these preliminary results with genetic analysis, starting to get an insight in population genetics in Nicaraguan honeybee colonies. Our limitations are due to our financial situation; however, we look forward to realize further research with approaches that are more detailed.

We are pleased to say that we have effected all necessary changes and corrections to your comments as you see in the following part of the responses to the reviewers and in the attached document named ‘Revised Manuscript with Track Changes’. As advised, we have also attached a separate document named ‘Manuscript’ that do not have track changes. Please do not hesitate to get back to us if there is any need to do so.

Best regards

The authors

Editor and Reviewer comments:

Journal Requirements:

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming.

Resp: Corrected. We ensured that our manuscript meets PLOS ONE's style requirements, including those for file naming.

2. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

Resp: Corrected. We provided all data and additional material as files of the supporting information (file name: “S1_Data.XLSX” and “S1_Fig.TIFF”).

3. We note that Figure 2 in your submission contain copyrighted images. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution.

Resp: Corrected. Dear editor, the photograph in figure 2 is not copyrighted, we took photographs of the samples in the laboratory, when we started the diagnostic procedures. When the article is published, we agree that our material can be used under Creative Commons Attribution License (CC BY 4.0).

This is the reference for figure 2 (Fig2.TIF):

Düttmann, C. (Photographer). 2013, November, León, Nicaragua. Honeybees sampled from one beehive in Boaco, sample no. A8-04-10-13. Laboratory CEVEDI – UNAN - León.

In the same way the new figure we provide for Supporting Information, as suggested by reviewer #1, is another photograph that is not copyrighted. We made the photograph using a stereoscope to determine the forewing length of the sampled honeybees and now we edited it to indicate measurement practice.

This is the reference for figure S1 (S1_Fig.TIFF):

Düttmann, C. (Photographer). 2016, September, León, Nicaragua. Right forewing of Apis mellifera detached for length measurement, Laboratory CEVEDI – UNAN - León.

4. We note that Figure 5 in your submission contain [map/satellite] images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

Resp: Corrected. We replaced the copyrighted map in figure 5 with a map from the following source: naturalearthdata.com. "All versions of Natural Earth raster + vector map data found on this website are in the public domain. You may use the maps in any manner, including modifying the content and design, electronic dissemination, and offset printing. The primary authors, Tom Patterson and Nathaniel Vaughn Kelso, and all other contributors renounce all financial claim to the maps and invites you to use them for personal, educational, and commercial purposes. No permission is needed to use Natural Earth. Crediting the authors is unnecessary.”

Comments to the Author:

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

5. Review Comments to the Author

Reviewer #1:

#General comments:

The manuscript discusses the Africanization level of honeybee (Apis mellifera 1 L.) from Nicaragua using phenotypic, molecular, behavior, and morphometric analyses. I believe that this manuscript can be published in this Journal, but it is necessary to do a major reviewer. The authors need to organize better all methods that they used. They also need to work more in the discussion.

Resp: Dear reviewer. Thank you for your consideration and your constructive suggestions. As you will see in the following points, we added the lacking information in method, we improved the reporting of the results and we enriched the discussion.

I suggest the authors add geometric morphometrics analysis to study shape and centroid size of the forewings because it can provide more accurate results than the traditional morphometrics that they used.

Resp: Dear reviewer, thank you for this suggestion. Of course, the study would have shown more results if we had carried out a geometric-morphometric analysis, but the main aim of the study was not to define subspecies by morphometric features, but to determine the degree of Africanization of the honeybees in Nicaragua by the length of right forewing as described by Rinderer et al. 1986. Phylogenetic lineages were determined by mtDNA analysis to confirm the results and to provide initial insights into honeybee populations.

#Abstract

Include the question and the main objective of the study in the first paragraph.

Resp: Corrected. We included the principle question and the main objective of the study in the first paragraph of the abstract.

Line 28: the authors used four methods to study honeybee populational variation. Please add the phenotypic and behavior in this line.

Resp: Corrected. We added that we related beehive characteristics and management to the results of morphometric and genetic analyses.

Line 28: Replace “approach” with approaches.

Resp: Corrected. We revised the manuscript and we corrected approach to approaches where it was necessary.

# Introduction

Line 51: replace “analysis” with analyses.

Resp: Corrected. We revised the manuscript and we corrected analysis to analyses where it was necessary.

Line 51-56: please extend this literature reviewer.

Resp: Corrected. We have extended our statement about the various morphometric and genetic approaches for the determination of Apis mellifera species, subspecies and populations with brief explanations and added further references to support them.

Line 102: Choose “behaviour” (UK) or “behavior” (USA) in your manuscript.

Resp: Corrected.

Line 115-118: which previous studies did you use to elaborate your hypotheses? Please, include them.

Resp: Corrected. We cited three studies from Costa Rica, Mexico and USA where the factors are discussed that support our hypothesis.

#Methods

Line 133: are bees collected in the same year?

Resp: Corrected. In the text, we included the period sampling was conducted, but the years are not relevant for statistical analysis.

Lines 145-154: I suggest this paragraph be a subtitle about behavior

Resp: Corrected. Subtitle added “Data collection about defensive behaviour and beehive management”

Lines 151-153: “Further information was recorded about the applied methods of multiplication of the beehives, queen bee replacement and the origin of the queen in the beehive”. What did you use to get this information?

Resp: Corrected. During the sampling, additional data were recollected using a questionnaire applied to the beekeepers who manage the selected apiary. The questionnaire was previously pilot-tested.

Line 163: Are the “146 samples” the beehives? Please, specify it.

Resp: Corrected. We improved the sentence for better understanding:

As the first feature in this study, phenotypic diversity of the bee colonies in each of the 146 hives was examined using photographic identification of the collected bees to document uniformity or diversity of the abdominal colour pattern.

Line 166-172: Please provide an image of the forewing indicating the measurements. You also have two variables: colour and wing length. How many bees of the same colour did you use to measure the forewings? It is interesting to investigate if there are differences among bee forewings with different behaviors.

Resp: Corrected. A) We added an image of the forewing measurement and included it to Supporting Information (S1_Fig.TIFF).

B) As Ruttner (1988) described, differences in colour among bee colonies is a bimodal variation, therefore our selection criteria to measure the forewing length, was based on the phenotypic diversity in the sample, including the same number of bees from each colour. We added information that is more detailed in methods.

C) We studied the association between the result of the mean forewing length of a bee colony and the beehive characteristics such as the origin of the hive and the defensive behaviour, as well as beekeeping practices. “There was no significant difference between the characteristics of the beehives and management activities according to the forewing length (ANOVA, origin of the colony: F3,146 = 0.13, p > 0.05; abdominal colour pattern: F1,146 = 0.55, p > 0.05; defensiveness of the beehive: F2,146 = 1.52, p > 0.05; queen replacement: F 1,146 = 0.48, p > 0.05; origin of the queen: F2,146 = 0.66, p > 0.05) (Fig 4)”

Line 172: Add the references of the software.

Resp: Corrected.

Line 176-182: why did you include samples in your work just for molecular analyses?

Resp: Dear reviewer, we included to our own sample collection during the study, samples from five different sentinel apiaries, which are of interest because of their geographic position “northern border to Honduras (1), southern border to Costa Rica (1), Island in the Nicaraguan Lake (1) and the South Caribbean Coast (2)”. These five locations could have been a possible provenience for different genetic findings. The samples were obtained directly from the Institute of Agricultural Protection and Health (IPSA), but included only few specimen for diagnostic procedure, therefore additional information was not available, such as phenotypic characters or forewing length measurement and on the other hand the lack of data from the questionnaire. However, the genetic data is important to get an insight into the existing lineages of Apis mellifera in Nicaragua.

Lines 183-186: why did you include a new hypothesis here? Which previous studies did you use to elaborate this hypothesis?

Resp: Dear reviewer, this is not a new hypothesis, here we underline our assumption that exchanging frames with brood between beehives to strengthen the weak hives could be a reason to find different mitotypes in the same beehive.

As Evans et al. indicated in “Standard methods for molecular research in Apis mellifera”, that “Due to the risk of drifting between colonies” (even if they mention natural drifting) “it is arguably worth sampling more than one individual to avoid mistakes in assigning colony heritage.” The impact of hive management on genetic pattern is still very important to investigate and the exchanging of brood combs is a widespread practice in Nicaragua.

#Results

Line 263: you also studied the abdominal colour

Resp: Dear reviewer, this is correct and we pointed it out in line 266.

Lines 274-283: the authors need to describe this procedure in the methods. How did you get this result?

Resp: Corrected. We added in methods “Management information was recorded from the beekeepers about the applied methods of multiplication of the beehives, queen bee replacement and the origin of the queen in the beehive. To collect data, the investigators developed a questionnaire that was pilot – tested before application.”

Line 294: Write the species name in italic. Please, look at all manuscript.

Resp: Corrected.

Lines 300-306: present only your results and move this information to discussion.

Resp: Corrected.

# Discussion

The authors need to discuss all result topics including the evolutionary history of Apis mellifera.

Resp: Dear reviewer, thank you for your comment. In the discussion section, we added a paragraph discussing the characteristics of beehives and the management of apiaries. In regards to the evolutionary history, we believe we have addressed all aspects relevant to our results in the discussion.

Lines 480-482: the authors agreed that geometric morphometrics can “offer more information about the evolutionary processes”. However, they used traditional morphometrics.

Resp: Dear reviewer, here we mentioned that for subspecies discrimination, it is necessary to extent the approach to geometric morphometrics, but our aim was to distinguish Africanized from European honeybee colonies. In this context (line 480-482), the discussion was about the longer forewings we found in bees with an A4 mitotype compared to those bees with A1 mitotype. If we consider that, the A1 mitotype that we sent to GenBank was identified as A. mellifera iberiensis and if we want to relate this finding to a morphometric approach, it would be necessary to apply geometric morphometrics to discriminate subspecies. Nevertheless, this was not the aim of our study; however, it is a relevant indication for future studies.

Reviewer #2:

The process of the Africanization of Apis mellifera and its occupation in the Americas is a fascinating phenomenon that still needs to be studied. Information on the genetic and morphological diversity of the honeybee in Nicaragua contributes to a greater understanding of this process.

Resp: Dear reviewer. Thank you for your comment; we share the same enthusiasm about Apis mellifera research and we want to fill an information gap in Central America with our first study.

Line 55: it would be interesting to mention what these recent techniques are about. Unfortunately, in the text, it was not clear.

Resp: Corrected. The same observation made reviewer #1. As you can see in our response to him/her, we have extended our statement about the various morphometric and genetic approaches for the determination of Apis mellifera species, subspecies and populations with brief explanations and added further references to support them.

Line 197: describe the modifications in the protocol.

Resp: Corrected. We added that we changed in the elution volume of the samples (from 200 µl to 50 µl for each sample), a common practice to decrease elution volume in genetic approaches. We decided to do this modification, because in our experience we noted better results with this elution volume. Decreasing elution volumes does not significantly compromise DNA yields (El-Mogy et al. The Effect of Elution Volume on DNA Quantity and Quality Using Norgen’s Saliva DNA Isolation Kit. 2016).

Lines 262-266: would fit better at the end of the text as a conclusion to the study.

Resp: Dear reviewer, thank you for this suggestion, but we regard these 5 lines as an introduction to the following results, a kind of summary which results are to be expected. However, we improved the conclusions concerning your observation.

General observations:

- Among a series of measures proposed by Ruttner in traditional morphometry to distinguish Apis mellifera lineages, do you think choosing only one linear measurement (affected by environmental factors) is enough to discriminate the groups? For example, according to Diniz-Filho et al., 2000 “Phylogenetically, the wing venation is more informative compared to the more environment–sensitive character categories of size, colour or pilosity”.

Resp: Dear reviewer, As we mentioned in our manuscript, now Line 112 - 114: “Even assuming that this technique (note: measurement of forewing length) provides only preliminary results [13] and that it is unsuitable for identifying subspecies, it does indicate the probability that a colony is Africanized [21].”

As we mentioned in our initial letter and in the response to reviewer #1, the main aim of our study was to identify the degree of Africanization in Nicaragua and the mitotype lineages distributed in the different zones where apiculture is most frequent.

Ruttner established an extended Morphometric Bee Data Bank in Oberursel, Germany (the Ruttner collection with international participation) to discriminate groups, we are far from this, but with our study we can contribute to future research hypothesis. Another factor is, that due to our financial and technical limitations, we chose the simplified technique described by Rinderer et al., using a single morphometric character to differentiate just between Africanized and European bee colonies. The principal result of the study Rinderer et al. “Field and Simplified techniques for identifying Africanized and European honeybees, 1986”, showed that the forewing length was “the best single character to distinguish between Africanized and European honeybees: 86% were correctly identified and no sample was misidentified; the remaining were unidentified.” Based on the morphometric results, we provided data that are more accurate for identifying lineages by mtDNA analysis and sequencing.

- In the text, the average length of the wing of European honeybees was not mentioned. Are they larger or smaller than the measurements of Africanized bees?

Resp: Corrected. We added that the overall average of the forewing length in European honeybees is �  = 9.20 mm, larger than the Africanized bees. Reference:

Rinderer TE, Sylvester HA, Brown MA, Villa JD, Pesante D, Collins AM, et al. Field and simplified techniques for identifying africanized and european honey bees. Apidologie. 1986;17: 33–48. doi:10.1051/apido:19860104

- Why was the entire sample of phenotypic variation not included in the molecular analyses?

Resp: Dear reviewer, all the samples with complete information were included and analysed, inclusive the phenotypic variation. However, additionally we included 14 honeybees from five sentinel apiaries provided by IPSA, the Institute of Agricultural Protection and Health (see answer to line 176 – 182 reviewer #1).

Despite the description of the criteria “...phenotypic diversity within the colony and mean length of the forewing.”, which factor determined this choice was unclear.

Resp: Corrected. After a carefully revision, we saw were the problem of understanding of our sample size for molecular diagnosis is. We changed the order of the paragraphs and we explained the procedure better.

For molecular diagnosis, honeybees from the saved 146 samples were selected from all inspected apiaries in study. The criteria heterogeneity of the abdominal colour within the colony and mean length of forewing were indicative for selection. However, first we had to mention that we decided to make PCR diagnosis with two worker bees instead of one (see Evans et al. Standard methods for molecular research in Apis mellifera). In our previous submitted manuscript, we explained this part in the second paragraph, so this led to confusion. The new redaction avoid this misunderstanding.

- Neighbor-joining is a clustering algorithm that clusters haplotypes based on genetic distance. Maximum likelihood can apply a model of sequence evolution and is better for building a phylogeny using sequence data. Neighbor-joining isn’t a phylogenetic analysis

Resp: Corrected. Thank you for this suggestion; we revised literature and we agree that ML method is accurate to determine evolutionary divergences. We substituted Neighbor-joining by the Maximum likelihood method and we updated section methods as well as the results for the phylogenetic tree.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

12 Apr 2022

Africanized honeybee population (Apis mellifera L.) in Nicaragua: Forewing length and mitotype lineages.

PONE-D-21-31700R1

Dear Dr. Flores,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Tiago M. Francoy, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

After carefully reading the new version of the manuscript and the anwers provided by the authors, I think all the requirements asked by the reviewers are satisfactory. It is my opinion that the manuscript is ready for acceptance and publication.

Reviewers' comments:

14 Apr 2022

PONE-D-21-31700R1

Africanized honeybee population (Apis mellifera L.) in Nicaragua: Forewing length and mitotype lineages.

Dear Dr. Flores:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Tiago M. Francoy

Academic Editor

PLOS ONE