*4.3. Data Processing, Statistical Analysis, and Reporting*

The known copy number of the standard for each virus, with 10-fold dilutions from 10−<sup>3</sup> to 10−7, were prepared and added in each RT-qPCR run. The exact number of RNA viral molecules in individual sample was calculated for positive samples from the standard curve for each of the four honeybee viruses.

The results for each sample were analyzed using MxPro-Mx3005P v4.10 software (Stratagene, La Jolla, CA, USA) and the exact copy number was determined from the standard curve. Results were expressed as number of detected viral copies in 5 μL of extracted RNA.

For sequencing purpose, the BQCV- and DWV-positive samples were amplified by using a specific method of reverse transcription and polymerase chain reaction (RT-PCR), as previously described [33]. Results were evaluated based on the size of RT-PCR products in the agarose gel as positive in the case of the expected product size: for BQCV, 770 nt, and for DWV, 504 nt [33]. In the case of a positive result, the selected RT-PCR products of a single virus were directly sequenced with the Sanger sequencing protocol, using the same primers as used for specific RT-PCR as described previously [10]. Individual sequences were analyzed using the DNASTAR 5.05 (Lasergen, WI, USA) program, and 26 positive samples of two viruses (BQCV *n* = 24 and DWV *n* = 2) were detected in bumblebees together with closely related sequences from GenBank and interpreted according to the results of the nucleotide sequence matching between honeybee and bumblebee samples. Multiple alignments were created using the program MEGA 6.06. Genetic distances were calculated from the alignment based on the Tamura three-parameter model, and phylogenetic trees were generated using the maximum likelihood (ML) statistical method implemented with the Tamura 3-parameter model with Gamma distribution [34]. The test of phylogeny was performed through 1000 bootstrap replicates. Only bootstrap values higher than 70% were presented on phylogenetic trees. The comparative analyses of collected samples at 27 locations and comparison with previously detected viruses in honeybees and with the most closely related sequences available in GenBank were performed.

#### **5. Conclusions**

Results of this research confirm that several honeybee viruses (ABPV, CBPV, BQCV, DWV) were found in wild bumblebees in Croatia. It can also be concluded that the presence of the examined viruses was higher in continental parts of the country compared with the Adriatic coast and islands. However, because the virus's presence was not studied in internal bumblebee tissues of separate bumblebee species, in further studies, samples of different bumblebee species will be tested individually to define species-specific prevalence and to evaluate the possible impact that those viruses may have on their bumblebee hosts. Moreover, due to the neighborhood and a good beekeeper connection via trade, virus genetic matches were determined in Croatia and Slovenia.

**Author Contributions:** Conceptualization, I.T.G. and I.T.; investigation, I.T.G., I.T. and L.Š.; writing original draft preparation, I.T.G.; writing—review and editing, I.T.G., I.T. and L.Š.; visualization, I.T. and I.T.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are presented in the manuscript. Twenty-six determined sequences from this study were deposited in GenBank with accession numbers MW488237-MW488262.

**Acknowledgments:** Authors thank the students of the Faculty of Veterinary Medicine University of Zagreb for their help with the collection and transportation of wild bumblebee samples.

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