Whole genome sequencing of Avian metapneumovirus type B genomes directly from clinical samples collected from chickens in live bird markets using multiplex tiling RT-PCR method

IntroductionAvian metapneumovirus (AMPV), a member of the family Pneumoviridae, genus Metapneumovirus possesses a non-segmented negative-sense RNA genome of approximately 13kb~15kb with eight genes (3’-N-P-M-F-M2-SH-G-L-5’) (1-3). AMPV causes rhinotracheitis in turkeys and swollen head syndrome (SHS) in chickens, mostly contributing to secondary bacterial infections leading to more severe symptoms in chickens (4). Different types of AMPVs have been classified based on the nucleotide sequence divergence of the attachment glycoprotein (G) and antigenic differences between strains (5). Types A and B are found all around the world, and type C was reported in North America, China, and South Korea, and in a retrospective study in France in 1990 (6). Type D was reported only once in a turkey flock in France in 1985 (5, 7). AMPV types A, B, and C have been detected in South Korea; the type A and B in poultry farms and the type C from pheasants in live bird markets (6, 8, 9). The study of AMPV infection is particularly difficult due to the transient nature of viral shedding in the host before the symptoms develop. In most cases, multiple blind passages of AMPV from the clinical samples are required for isolation and identification (4, 5). Therefore, a limited number of AMPV sequences are available in the NCBI GenBank database. As of November 29, 2022, there are only three complete genome sequences of AMPV type B (LN16, VCO3/60616, and Hungary/657/4; MH745147, AB548428, and MN729604 respectively) (3, 10, 11). The tiling amplicon method has been proven to be efficient and prolific in producing whole genome sequences directly from clinical samples during the Covid-19 pandemic (12). In a previous study, the multiplex tiling RT-PCR method was applied to enrich the genetic material of Zika virus and Rabies lyssavirus directly from brain tissue samples, and it enables to sequence of the full genome of low titer samples containing as low as 50 genome copies in a reaction (13, 14).In this study, we detected 6 AMPVs from chickens in live bird markets (LBM) in South Korea during 2019-2022. For whole genome sequencing of AMPVs, we developed a PCR primer panel to efficiently amplify the complete coding region of AMPV type B using the multiplex tiling RT-PCR method. We successfully obtained the complete coding region of AMPVs using Illumina next-generation sequencing (NGS) and conducted comparative a phylogenetic analysis to analyze the genetic relatedness of AMPVs from LBMs in Korea with other AMPVs.Material and MethodsSample collectionA total of 138 slaughtered chickens were purchased from poultry meat vendors in LBMs in Korea during the period of 2019 to 2022. The nasal turbinate or whole beaks of the chickens were collected and washed with phosphate-buffered saline (PBS) using enough to fully immerse the samples. Total RNA was extracted from the nasal turbinate wash samples using the Qiagen RNeasy mini kit (Hilden, Germany) and used to detect both avian metapneumovirus type A and B using the real-time qRT-PCR as previously described (8). All positive nasal turbinate wash samples were inoculated to Vero cells for virus isolation.Primer design and tiling amplicon PCRTo design a tiling amplicon PCR primer panel, full genome sequences of AMPV type B available in the NCBI GenBank (LN16, VCO3/60616, and Hungary/657/4; Accession no. MH745147, AB548428, and MN729604 respectively) were downloaded and aligned using the MAFFT 1.4.0 program on Geneious Prime software (https://www.geneious.com) (3, 10, 11). Primers were designed to amplify 380-420bp region with about 100bp overlap, and 43 sets of designed primers were pooled into two pools according to the primer design output by Primal Scheme (http://primalscheme.com) (Table S2) (13). Three of six positive samples with lower cycle threshold (Ct) values (ranging from 25.77 to 30.6) were selected for subsequent multiplex RT-PCR (Table 1). cDNA was synthesized from the RNA using the LunaScript RT SuperMix kit (NEB, Massachusetts, United States) following the manufacturer’s instructions. For the multiplex PCR assay, an equal volume of each 100uM primer stock was pooled together as designated as pool 1 and 2 (Table S2). The PCR mixture was prepared by mixing 12.5ul Q5 Hotstart 2X Master Mix (NEB, Massachusetts, United states), 0.015uM of each primer pool, 5ul of template cDNA, and nuclease-free water up to 25ul. The reaction mixture was prepared for each pool as it induces consistency between reactions and PCR assay was performed. PCR amplification conditions were: 98℃ for 15s, followed by 35 cycles of 95℃ for 15s and 63℃ for 5min (12). The PCR products were then visualized by electrophoresis on 2% agarose gels showing around 400bp amplicons (13). PCR products of both primer pools were combined and used for library preparation after purification using the Qiaquick PCR purification kit (Qiagen, Hilden, Germany). Library was prepared using the TruSeq Nano DNA kit (Illumina, California, United States) to produce 2x151 bp paired-end reads. Library preparation and sequencing on NextSeq 500 sequencing system (Illumina, California, United States) was done by LAS (Gimpo, Republic of Korea). Assembly and phylogenetic analysisRaw reads were trimmed of adapters and low-quality bases using BBDuk version 38.84 by setting the minimum quality to 20 (15). De novo and directed assemblies of genome sequences were performed. Trimmed reads were mapped to the LN16 virus genome (GenBank accession number: MH745147) using the Minimap 2.24 (https://github.com/lh3/minimap2) with default options and visualized on Geneious Prime software. Trimmed reads were assembled de novo using the SPAdes assembler 3.15.5. The assembled genome sequences produced by reference guided genome assembly and de novo assembly were combined to generate the final consensus genome sequences. A total of 48 glycoprotein (G) sequences of AMPV type B were downloaded from the GenBank database and aligned using the MAFFT Multiple Sequence Alignment software v7.450 for phylogenetic analysis (16). A maximum-likelihood (ML) phylogenetic tree was reconstructed using the RAxML GUI 2.0 (https://antonellilab.github.io/raxmlGUI/) with a rapid bootstrap option set to 1,000 (17).Descriptive resultsA total of 6 out of 138 (4.34%) chicken nasal turbinate wash samples tested positive for AMPV type B by real-time qRT-PCR (Table S1 and Supplementary Figure 1). The inoculated Vero cells did not exhibit any cytopathic effect (CPE), suggesting it may require additional passages for virus isolation (data not shown). We successfully obtained complete coding genome sequences (CDS) of all three nasal turbinate samples (AMPV/B/Korea/N19-29/2019, AMPV/B/Korea/N21-83/2021, and AMPV/B/Korea/N21-41/2021; hereafter N19-29, N21-83, and N21-41, respectively) (Table 1). Our tiling amplicon PCR reactions showed a fairly small amplification bias when clinical samples with relatively low Ct values, 25.77 and 28.67, were used. We were able to assemble complete CDS of N19-29 and N21-41 viruses using the tiling amplicon PCR primer panel designed in this study coupled with illumina NGS. However, the N21-83 virus which showed a higher Ct value (30.6) had two short gaps in the initial genome assembly results, including 5bp in F gene and 100bp in L gene. The gaps were covered with two sanger sequencing reads using the same primers used in this study. Since the two fragments were amplified using the same primers, the outcome of low amplification efficiency was not likely due to primer-template mismatch. We assume that competitive inhibition between primers may have decreased PCR efficiency of the fragments. Phylogenetic analysis of G gene sequences revealed that the viruses sequenced in this study belong to the AMPV type B and showed the closest genetic relationship with the SC1509 (GenBank Accession no. DI187010.1) virus (sequence identity: 95.74-96.14%) from Korea in 2009 and the LN16 virus (sequence identity: 95.42-95.82%) identified from China in 2016 (Figure 1). The N19-29, N21-83, and N21-41 viruses did not cluster with the commercial live attenuated AMPV type B vaccine strains in phylogenies (Figure 1). In conclusion, we developed the tiling amplicon PCR method for genome sequencing and successfully sequenced three AMPVs directly from clinical samples. The multiplex tiling RT-PCR and NGS approach developed in this study has the potential to be implemented in a diagnostic setting, providing a rapid and reliable method for complete genome sequencing and molecular epidemiological study of AMPV from clinical samples. In addition, the complete CDS of AMPVs established in this study would be useful as reference data for future investigations on AMPVs. The relative ease of acquiring complete CDS directly from clinical samples without a labor-intensive adaptation process in cell culture will help further the diversity of the AMPV genome database. In addition, since the live bird markets have been recognized as a reservoir, amplifier, and source of avian viruses (18), genomic surveillance of AMPVs in LBMs should be enhanced for monitoring of further evolution and spread of the AMPVs.