Tóm tắt
Công nghệ giải trình tự thế hệ mới (NGS) đã cách mạng hóa phương pháp phân tích gen 16S rRNA trong các nghiên cứu sinh thái vi sinh vật. Hiện có ba công nghệ chủ đạo cho việc giải trình tự 16S quy mô lớn, mỗi công nghệ đều bị ảnh hưởng bởi các yếu tố như sai sót giải trình tự, tính đặc hiệu của mồi khuếch đại và độ dài đoạn đọc, vốn có thể làm sai lệch bức tranh thực tế của quần xã vi sinh vật. Nghiên cứu này so sánh hiệu quả của việc sử dụng đoạn đọc ngắn (vùng V1-V3) và đoạn đọc dài gần toàn bộ gen (vùng V1-V8) của gen 16S rRNA trong việc phân tích quần xã vi sinh vật dạ cỏ bò, một hệ sinh thái có độ đa dạng cao, nhằm đánh giá tác động của lựa chọn công nghệ lên hồ sơ phân loại phát sinh loài. Trình tự vùng V1-V3 được thực hiện bằng phương pháp đọc ngắn ghép cặp trên nền tảng Illumina MiSeq, trong khi trình tự vùng V1-V8 được tạo ra bằng công nghệ đọc dài “circular consensus” trên thiết bị Pacific Biosciences RSII. Kết quả cho thấy cả hai nền tảng đều mang lại các chỉ số tương đồng về số lượng đơn vị phân loại hoạt động (OTU), độ phong phú loài, chỉ số bao phủ Good (Good’s coverage) và chỉ số đa dạng Shannon. Tuy nhiên, phân tích vùng V1-V8 cho thấy sự gia tăng đáng kể về mức độ phong phú tương đối của một số bậc phân loại, đặc biệt là ngành Proteobacteria và Verrucomicrobia (P < 0,05). Độ chính xác trong phân loại định danh cũng cao hơn khi sử dụng đoạn đọc dài. Phân tích ma trận khoảng cách UniFrac bằng phương pháp phân cụm UPGMA có kiểm định jackknife cũng cho thấy sự khác biệt rõ rệt giữa các quần xã được xác định bởi hai phương pháp. Các kết quả này củng cố quan điểm rằng đoạn đọc dài hơn cung cấp độ phân giải phát sinh loài vượt trội, điều mà các đoạn gen 16S rRNA ngắn hơn có thể không đạt được. Nghiên cứu này trên quần xã vi khuẩn dạ cỏ bò cho thấy việc sử dụng đoạn đọc 16S gần toàn bộ chiều dài là cần thiết cho các phân tích chuyên sâu, hoặc để xây dựng cơ sở dữ liệu tham chiếu đặc thù cho từng môi trường. Cơ sở dữ liệu này sau đó có thể được ứng dụng để cải thiện độ chính xác của các phân tích dựa trên công nghệ đọc ngắn khi chi phí là một rào cản đối với việc giải trình tự 16S toàn bộ gen.
1. Giới thiệu
Công nghệ giải trình tự thế hệ mới (NGS) đã trở thành một công cụ trung tâm trong các nghiên cứu sinh thái vi sinh vật, đặc biệt là các phương pháp không dựa trên nuôi cấy, vốn dựa trên phân tích phát sinh loài phân tử của gen RNA ribosome tiểu đơn vị nhỏ (gen 16S rRNA). Khả năng liên kết các xu hướng ở cấp độ loài hoặc chi với các thông số vật chủ/môi trường thông qua việc lập hồ sơ 16S đã được chứng minh là một công cụ mạnh mẽ. Các nghiên cứu dựa trên khuếch đại gen này phụ thuộc vào khả năng bắt cặp của mồi khuếch đại với các vùng bảo tồn bao quanh các vùng siêu biến (V1-V9) của gen 16S rRNA. Tuy nhiên, các sai sót cố hữu trong quá trình giải trình tự amplicon và việc lựa chọn vùng siêu biến có thể làm sai lệch các đánh giá định lượng về động thái của quần xã vi khuẩn, ngay cả khi đạt được độ sâu giải trình tự cao. Vấn đề này đã và đang là một trọng tâm nghiên cứu, với mục tiêu giảm thiểu những hạn chế của việc ước tính đa dạng vi sinh vật dựa trên PCR. Việc lựa chọn mồi cũng ảnh hưởng đến hồ sơ phát sinh loài biểu kiến, và hiện vẫn chưa có sự đồng thuận về việc nên sử dụng một hay nhiều vùng siêu biến nào của gen 16S rRNA cho phân tích quần xã. Theo tiêu chuẩn vàng, gen 16S rRNA toàn bộ chiều dài (khoảng 1500 bp) có thể được sử dụng để định danh phân loại chính xác, mặc dù hạn chế về chi phí thường dẫn đến việc chỉ nhắm mục tiêu một hoặc một vài vùng siêu biến. Điều này cho thấy rõ ràng rằng việc lựa chọn vùng siêu biến có thể ảnh hưởng đến kết quả phân loại, độ phong phú và tính đa dạng của OTU. Các nghiên cứu này làm sáng tỏ rằng việc sử dụng trình tự toàn bộ chiều dài, hoặc tiến hành các nghiên cứu sơ bộ để xác định vùng siêu biến tối ưu cho môi trường sinh thái đang được nghiên cứu, sẽ cải thiện hiệu quả của kỹ thuật lập hồ sơ dựa trên 16S. Các đoạn đọc gần toàn bộ chiều dài chất lượng cao có thể được tạo ra trên thiết bị Pacific Biosciences RSII, sử dụng công nghệ giải trình tự SMRT đọc dài ở chế độ “circular consensus”. Tuy nhiên, chi phí cho mỗi đoạn đọc cao hơn đáng kể so với các nền tảng đọc ngắn, điều này có thể hạn chế việc ứng dụng công nghệ này cho các so sánh quy mô lớn trên nhiều mẫu. Do hạn chế về chi phí, các nghiên cứu hiện nay thường tập trung vào các vùng ngắn hơn, như V1-V3 hoặc V3-V5, được giải trình tự trên nền tảng Illumina MiSeq hoặc Life Technologies Ion Torrent. Hệ sinh thái vi sinh vật của dạ cỏ gia súc đã được nghiên cứu rộng rãi và sở hữu cấu trúc phức tạp, đa dạng, cần thiết cho vật chủ trong quá trình tiêu hóa vật liệu thực vật. Các nghiên cứu gần đây đã tập trung vào việc mô tả đặc điểm hệ vi sinh vật dạ cỏ giữa các động vật khác nhau về hiệu quả sử dụng thức ăn. Những nghiên cứu này đã mô tả quần xã vi khuẩn bằng cách giải trình tự vùng V1-V3 hoặc các vùng V1/V3 riêng lẻ của gen 16S rRNA. Phân tích các vùng siêu biến này được báo cáo là có khả năng chứng minh các xu hướng, với một số thay đổi đáng kể ở các đơn vị phân loại cụ thể, nhưng nhìn chung, các nghiên cứu này còn thiếu độ phân giải phát sinh loài. Các nghiên cứu trước đây của chúng tôi đã chứng minh mối liên hệ ở cấp độ chi và loài của vi khuẩn với sự khác biệt lớn về hiệu quả sử dụng thức ăn ở bò đực, nhưng cũng lưu ý rằng để xác định chính xác hơn các vi khuẩn liên quan đến hiệu quả sử dụng thức ăn, cần có các phân tích với độ phân giải cao hơn, chẳng hạn như phân tích metagenome (giải trình tự shotgun), xây dựng một cơ sở dữ liệu tham chiếu đặc thù cho môi trường nghiên cứu, hoặc giải trình tự toàn bộ gen 16S rRNA. Dữ liệu trước đây của chúng tôi về tác động của sự biến đổi hệ vi sinh vật đối với hiệu quả sử dụng thức ăn ở bò đực đã chỉ ra sự cần thiết phải tăng cường độ phân giải của phân tích. Do đó, nghiên cứu này nhằm mục đích đối chiếu hồ sơ vi sinh vật của quần xã vi khuẩn dạ cỏ được tạo ra từ vùng V1-V3 và V1-V8 của gen 16S rRNA, sử dụng nền tảng Illumina MiSeq (V1-V3) và thiết bị Pacific Biosciences RSII (V1-V8). Chúng tôi đặt giả thuyết rằng sẽ có sự khác biệt rõ rệt về khả năng định danh đơn vị phân loại, số lượng và độ phong phú OTU, kết quả phân tích phát sinh loài, và hiệu quả sử dụng cơ sở dữ liệu tham chiếu giữa phương pháp giải trình tự một phần và gần toàn bộ gen 16S. Chúng tôi cũng mong muốn xác định liệu một cơ sở dữ liệu các trình tự 16S rRNA gần toàn bộ chiều dài từ dạ cỏ gia súc, được tạo ra bằng cách giải trình tự vùng V1-V8, có hỗ trợ việc xác định quần xã vi sinh vật chính xác hơn so với các cơ sở dữ liệu rRNA công cộng hay không, khi được sử dụng làm tham chiếu cho phân loại phát sinh loài của vùng V1-V3. Cuối cùng, nghiên cứu cũng nhằm xác định số lượng trình tự cần thiết để đạt được độ bao phủ toàn diện đối với tất cả các đơn vị phân loại được phát hiện bởi mỗi phương pháp khuếch đại (đoạn dài và đoạn ngắn).
2. Vật liệu và Phương pháp
2.1 Thiết kế thí nghiệm và thu thập mẫu dạ cỏ
Nghiên cứu này đã được phê duyệt bởi Ủy ban Chăm sóc và Sử dụng Động vật của Trung tâm Nghiên cứu Động vật Lấy Thịt Hoa Kỳ (U.S. Meat Animal Research Center). Hiệu quả sử dụng thức ăn được xác định theo quy trình của Myer et al. (2015a) và được sử dụng làm tiêu chí lựa chọn. Ba con bò đực có kiểu hình hiệu quả sử dụng thức ăn tương đương và độ lệch giữa chúng là thấp nhất đã được chọn và lấy mẫu cho nghiên cứu này.
[Đánh dấu để chèn hình ảnh về thiết kế thí nghiệm hoặc quy trình lấy mẫu]
2.2 Chiết xuất DNA, khuếch đại và giải trình tự
DNA tổng số được chiết xuất từ các mẫu dịch dạ cỏ bằng phương pháp nghiền hạt lặp lại kết hợp với cột (RBB+C). Thư viện DNA cho vùng V1–V3 của gen 16S rRNA được chuẩn bị bằng phương pháp khuếch đại PCR, sử dụng cặp mồi phổ biến đã được điều chỉnh 27F và 519R, có tích hợp trình tự adapter TruSeq® và chỉ số (index) cho nền tảng Illumina. Việc khuếch đại PCR và chuẩn bị thư viện DNA của vùng V1–V8 được thực hiện bằng cặp mồi phổ biến 27F và 1392R cho thiết bị Pacific Biosciences. Quá trình khuếch đại PCR được thực hiện với 20 chu kỳ cho cả hai vùng gen, sử dụng nhiệt độ bắt cặp là 58°C. Sản phẩm PCR được tinh sạch và các thư viện được định lượng. Các thư viện amplicon sau đó được giải trình tự bằng kit 2×300 bp, v3 (600 chu kỳ) trên nền tảng Illumina MiSeq® (đối với vùng V1–V3) và trên thiết bị Pacific Biosciences RSII (đối với vùng V1–V8). DNA metagenomic được phân lập, cắt đoạn và sử dụng để tạo thư viện TruSeq® PCR Free cho giải trình tự trên nền tảng Illumina NextSeq 500®. Tổng cộng, 147.995.952 đoạn đọc shotgun (tương đương 22,2 Gb dữ liệu) đã được tạo ra từ một lần giải trình tự.
2.3 Xử lý và phân tích trình tự đọc
Toàn bộ dữ liệu trình tự được xử lý bằng bộ phần mềm QIIME-1.9.1 và Mothur phiên bản 1.36.1. Đối với dữ liệu từ Illumina, các đoạn đọc ghép cặp (paired-end reads) được nối lại và lọc dựa trên điểm chất lượng (≥Q30). Các trình tự có độ dài ngắn hơn 400 bp bị loại bỏ. Đối với dữ liệu từ Pacific Biosciences, điểm chất lượng bằng 0 được coi là vị trí base không xác định và các trình tự được lọc dựa trên điểm chất lượng (≥Q30). Các trình tự có độ dài ngắn hơn 1200 bp bị loại bỏ. Đối với tất cả các đoạn đọc, các đoạn homopolymer dài hơn 7 bp bị loại bỏ và các trình tự chimeric được sàng lọc. Tất cả các trình tự đã được làm sạch được phân loại định danh bằng cơ sở dữ liệu gen 16S rRNA Greengenes, phiên bản 13_8. Các đơn vị phân loại hoạt động (OTU) được xác định bằng thuật toán uclust với ngưỡng tương đồng 97% (độ không tương đồng 0,03). Các OTU chỉ xuất hiện một lần (singletons) bị loại khỏi phân tích. Dựa trên phân tích đường cong bão hòa (rarefaction curves), dữ liệu được chuẩn hóa bằng cách lấy ngẫu nhiên một tập con gồm 25.000 trình tự/mẫu (dữ liệu MiSeq) và 40.000 trình tự/mẫu (dữ liệu Pacific Biosciences) để đảm bảo tính so sánh. Một cây phát sinh loài được xây dựng để xác định các chỉ số đa dạng alpha và beta. Các đoạn đọc metagenomic shotgun cũng được làm sạch tương tự và ánh xạ (map) vào các trình tự 16S rRNA đồng thuận (consensus sequences) của vùng V1-V3 và V1-V8. Các đoạn đọc ánh xạ vào các vùng siêu biến tương ứng đã được sử dụng để phân tích và phân loại.
2.4 Phân tích thống kê
Sự khác biệt về các chỉ số đa dạng trung bình và độ chính xác phân loại giữa các nhóm (phân tích vùng V1-V3 so với V1-V8) được đánh giá bằng phân tích phương sai một yếu tố (ANOVA) và kiểm định HSD của Tukey. Sự khác biệt có ý nghĩa thống kê được xác định ở mức P < 0,05. Đa dạng beta được ước tính bằng cách xác định khoảng cách UniFrac có trọng số và không trọng số giữa các mẫu, sử dụng cây phát sinh loài đã xây dựng và phần mềm QIIME. Tính có ý nghĩa của các cây được xây dựng được kiểm định bằng UniFrac significance test và P-test. Các mẫu được phân cụm dựa trên ma trận khoảng cách giữa chúng bằng thuật toán UPGMA. Độ tin cậy của cây UPGMA được ước tính bằng phương pháp jackknife. Phân tích thành phần chính (PCoA) được thực hiện dựa trên ma trận khoảng cách và kết quả phân tích UniFrac có trọng số và không trọng số. Đường cong bão hòa (rarefaction curves) được tính toán.
3. Kết quả
3.1 Các chỉ số giải trình tự
Kết quả sơ bộ từ phân tích mẫu dạ cỏ bò cho thấy 25.000 đoạn đọc/mẫu là đủ để bao phủ quần xã vi sinh vật khi sử dụng vùng V1-V3 của gen 16S rRNA cho phân tích. Độ dốc của đường cong bão hòa (rarefaction curve), biểu thị mối quan hệ giữa số lượng trình tự và số lượng OTU quan sát được, đã được so sánh giữa hai bộ dữ liệu amplicon.
[Đánh dấu để chèn Hình 1: Đường cong Rarefaction OTU cho trình tự Illumina (V1-V3) và PacBio (V1-V8)]
Độ dốc tại điểm 25.000 đoạn đọc của đường cong V1-V3 được ngoại suy sang đường cong V1-V8, cho thấy khoảng 40.000 đoạn đọc/mẫu đối với vùng V1-V8 là đủ cho phân tích so sánh. Do đó, các mẫu đã được chuẩn hóa (rarefied) xuống còn 25.000 và 40.000 trình tự/mẫu tương ứng cho vùng V1-V3 và V1-V8, được xem là đại diện cho độ sâu giải trình tự đầy đủ và được sử dụng cho các phân tích so sánh. DNA tổng số được chiết xuất từ dịch dạ cỏ của ba con bò đực. Cả ba mẫu DNA đều được khuếch đại bằng cả hai bộ mồi V1-V3 và V1-V8 trong các phản ứng riêng biệt (tổng cộng 6 phản ứng; ba phản ứng với V1-V3 và ba phản ứng với V1-V8). Sản phẩm khuếch đại được xử lý thành thư viện cho giải trình tự trên MiSeq hoặc RSII tương ứng. Sau quá trình lọc chất lượng và kiểm tra tính đặc hiệu của trình tự 16S, thu được 890.086 đoạn đọc từ nền tảng MiSeq và 138.180 đoạn đọc từ nền tảng RSII. Đối với mỗi cá thể động vật, số lượng đoạn đọc này được chuẩn hóa thành 25.000 (MiSeq) hoặc 40.000 (RSII).
[Đánh dấu để chèn Bảng 1: Thống kê đa dạng giữa các đoạn đọc từ các nhóm mẫu (V1-V8 vs V1-V3)]
Số lượng OTU trung bình từ các nhóm đã chuẩn hóa không khác biệt đáng kể (P > 0,05) giữa hai bộ dữ liệu amplicon. Tương tự, không có sự khác biệt đáng kể nào được tìm thấy trong các chỉ số phong phú (P > 0,05) giữa số lượng OTU. Mặc dù có một số biến động nhất định, các chỉ số đa dạng của mỗi nhóm sau chuẩn hóa không cho thấy sự khác biệt có ý nghĩa thống kê. Mỗi nhóm đều được bao phủ đầy đủ, lên tới 98,85% trong nhóm V1-V8 đã chuẩn hóa và 98,47% trong nhóm V1-V3 đã chuẩn hóa.
3.2 Thành phần phân loại
Độ phong phú tương đối của các ngành và chi vi khuẩn chiếm ≥ 1% tổng số đoạn đọc được trình bày trong Hình 2. Dữ liệu chuẩn hóa từ MiSeq (V1-V3) được phân loại thành 20 ngành, 28 lớp, 39 bộ, 64 họ và 73 chi. Ngành Bacteroidetes (68,64%) và Firmicutes (21,58%) chiếm ưu thế trong tập dữ liệu ở cấp độ ngành trong mỗi nhóm đã chuẩn hóa. Các đoạn đọc V1-V8 đã chuẩn hóa được phân loại thành 18 ngành, 27 lớp, 44 bộ, 65 họ và 81 chi. Tương tự như kết quả từ đoạn khuếch đại ngắn hơn, ngành Bacteroidetes (62,18%) và Firmicutes (23,38%) cũng chiếm ưu thế trong tập dữ liệu này.
Trong dữ liệu từ đoạn khuếch đại 16S gần toàn bộ chiều dài, tỷ lệ dữ liệu không được đại diện bởi hai ngành chiếm ưu thế này thấp hơn đáng kể. Độ phong phú tương đối của ngành Proteobacteria khác biệt đáng kể (P<0,05) giữa đoạn khuếch đại V1-V8 và V1-V3 (9,59% so với 0,51%). Cũng có sự khác biệt đáng kể giữa các đoạn khuếch đại đối với ngành Verrucomicrobia (0,3% so với 0,03%). Ngoài ra, tỷ lệ các đoạn đọc không được gán tên (unassigned reads) lớn hơn nhiều trong dữ liệu trình tự V1-V3 (7% so với 3% ở V1-V8). Việc phân loại đến cấp độ chi khả thi đối với 77% trình tự từ đoạn khuếch đại ngắn hơn và 72% từ đoạn khuếch đại dài hơn. Gần như tất cả (98% hoặc 97%) các đơn vị phân loại thuộc ngành Bacteriodetes đều bao gồm chi Prevotella, cho thấy sự phổ biến tổng thể của chi này có thể lên tới 67% trong các mẫu dạ cỏ này. Chi phong phú thứ hai, dựa trên phân tích đoạn khuếch đại dài hơn, là các thành viên của họ Succinivibrionaceae, mặc dù phân tích đoạn khuếch đại ngắn hơn không xếp họ này vào nhóm phong phú nhất. Để đánh giá khả năng sai lệch do mồi (primer bias), các đoạn đọc V1-V3 nhân tạo (in silico V1-V3) đã được tạo ra bằng cách cắt ngắn các đoạn đọc V1-V8 từ PacBio, sau đó tiến hành phân loại. Hồ sơ phân loại thu được cho thấy sự khác biệt giữa các đoạn PacBio V1-V3 nhân tạo và MiSeq V1-V3, cũng như với đoạn khuếch đại PacBio V1-V8. Kết quả cho thấy dấu hiệu của sai lệch do mồi, khi hồ sơ phân loại cấp ngành của các đoạn PacBio V1-V8 và PacBio V1-V3 nhân tạo có sự tương đồng cao về thành phần, nhưng khác biệt so với đoạn khuếch đại MiSeq V1-V3, đặc biệt về độ phong phú và khả năng phân loại của ngành Proteobacteria. Để kiểm tra sâu hơn khả năng sai lệch do mồi, giải trình tự trực tiếp DNA metagenomic đã được sử dụng để loại bỏ các sai lệch liên quan đến mồi. Trong dữ liệu này, DNA metagenomic được ánh xạ vào các vùng 16S rRNA V1-V3 và V1-V8 đồng thuận. Độ phong phú tương đối của các ngành và chi vi khuẩn chiếm ≥ 1% tổng số đoạn đọc được trình bày trong Hình 3.
Ở cấp độ ngành, mặc dù có sự khác biệt đáng chú ý giữa hồ sơ metagenomic và hồ sơ từ đoạn khuếch đại, đặc biệt là khả năng ước tính quá mức ngành Bacteroidetes từ vùng MiSeq V1-V3, dữ liệu từ đoạn khuếch đại V1-V3 vẫn khác biệt rõ ràng so với các hồ sơ còn lại về độ phong phú của Proteobacteria. Số lượng đoạn đọc không được gán tên đã giảm trong bộ dữ liệu 16S metagenomic. Các mối tương quan tương tự cũng được thể hiện ở cấp độ chi, tái khẳng định sự khác biệt về độ phong phú tương đối ở cấp độ chi. Sự khác biệt quan sát được về độ phong phú của các đoạn đọc không được phân loại giữa các vùng siêu biến kết hợp V1-V3 và V1-V8 đã được nghiên cứu sâu hơn bằng cách so sánh tỷ lệ các trình tự 16S rRNA do PacBio và Illumina tạo ra được phân loại và gán tên thành công ở năm cấp độ phân loại. Phân tích này cho thấy tỷ lệ phân loại thành công của các vùng siêu biến kết hợp V1–V8 lớn hơn đáng kể (P < 0,05) so với các vùng siêu biến kết hợp V1–V3, ở tất cả năm cấp độ phân loại, khi sử dụng Cơ sở dữ liệu gen 16S rRNA Greengenes.
Đúng như dự đoán, đối với tất cả các vùng siêu biến được phân tích, tỷ lệ phân loại thành công giảm dần khi tiến gần đến cấp độ chi; trong đó, tỷ lệ thành công cao nhất được ghi nhận ở cấp độ ngành và thấp nhất ở cấp độ chi.
3.3 Đa dạng Beta
Phân tích đa dạng beta được sử dụng để đánh giá mức độ tương đồng hoặc khác biệt về cấu trúc quần xã vi sinh vật giữa các mẫu. Các chỉ số dựa trên thông tin phát sinh loài cho phép ghi nhận những thay đổi trong thành phần và mối quan hệ tiến hóa của quần xã, không chỉ đơn thuần là sự khác biệt về thành viên giữa các quần xã. Do đó, các mẫu dạ cỏ được so sánh với nhau thông qua ma trận khoảng cách UniFrac sử dụng phương pháp phân cụm UPGMA có kiểm định jackknife, dựa trên dữ liệu từ đoạn khuếch đại dài hơn hoặc ngắn hơn.
Mặc dù các nút trên cây phân loại có thể ít được hỗ trợ hơn (lower bootstrap/jackknife support), việc phân cụm các mẫu cho thấy sự khác biệt rõ ràng khi so sánh kết quả từ vùng 16S V1-V8 và V1-V3 (Hình 5) trong cả ma trận khoảng cách không trọng số và có trọng số. Chiều dài nhánh cũng thể hiện sự khác biệt trong thước đo UniFrac không trọng số. Hơn nữa, mặc dù mỗi bộ dữ liệu amplicon chỉ bao gồm ba mẫu, các mẫu tương ứng từ mỗi cá thể lại phân tách khác nhau trong phân tích UniFrac có trọng số, tùy thuộc vào đoạn khuếch đại được sử dụng để phân tích. Để làm nổi bật sự tương phản này, phân tích thành phần chính (PCoA) đã được thực hiện, nhằm minh họa rõ hơn sự khác biệt so với phân cụm UPGMA (Hình 6).
Các kiểm định ý nghĩa UniFrac không trọng số và có trọng số, cùng với P-test, đã được sử dụng để xác định sự khác biệt giữa các quần xã vi khuẩn bằng cách tính toán khoảng cách giữa các mẫu dạ cỏ. Kết quả cho thấy các quần xã vi sinh vật dạ cỏ được xác định bằng các vùng trình tự 16S gần toàn bộ chiều dài và một phần là khác biệt đáng kể (P < 0,01).
4. Thảo luận
Việc lựa chọn đoạn gen 16S rRNA để giải trình tự trong nghiên cứu hồ sơ quần xã vi sinh vật phụ thuộc vào nhiều yếu tố, bao gồm ngân sách, đặc điểm môi trường nghiên cứu, độ bao phủ mong muốn đối với các nhóm sinh vật, và nền tảng giải trình tự. Dữ liệu thu được từ mỗi hệ sinh thái cụ thể có thể bị ảnh hưởng bởi nhiễu và sai lệch do việc lựa chọn vùng siêu biến. Hơn nữa, mức độ đa dạng trình tự trong gen 16S rRNA thay đổi giữa các vùng, điều này cũng có thể ảnh hưởng đến việc lập hồ sơ chính xác quần xã vi khuẩn. Trọng tâm của nghiên cứu này là đánh giá hiệu quả của việc sử dụng trình tự 16S một phần so với gần toàn bộ gen trong phân tích quần xã vi khuẩn dạ cỏ bò, chứ không nhằm so sánh trực tiếp các công nghệ giải trình tự. Khi kiểm tra đa dạng alpha của các mẫu từ hai đoạn khuếch đại khác nhau, kết quả cho thấy không có sự khác biệt giữa bất kỳ chỉ số nào được đo lường, đặc biệt là số lượng và độ phong phú của OTU (Bảng 1). Những dữ liệu này cũng xác nhận tính hợp lý của phương pháp so sánh dựa trên đường cong bão hòa OTU. Ngoài ra, nghiên cứu có thể thu hồi khoảng 98% tất cả các OTU được tính toán ở ngưỡng không tương đồng 0,03, như được xác định bằng ước tính độ bao phủ của Good, điều này chứng minh độ sâu giải trình tự đầy đủ đã được thiết lập trong các nghiên cứu trước đây. Hồ sơ phân loại vi khuẩn tương đối được tạo ra bởi các nhóm đoạn khuếch đại phù hợp với các nghiên cứu trước đây về quần xã vi sinh vật dạ cỏ, nơi độ phong phú hiện tại của Prevotella được ước tính từ 42% đến 60%, dựa trên trình tự gen 16S rRNA của vi khuẩn. Cũng phù hợp với nghiên cứu hiện tại, Prevotella đã được báo cáo là chi dạ cỏ phong phú nhất. Điều này quan trọng trong việc xác minh độ tin cậy của bộ dữ liệu, đặc biệt liên quan đến vùng V1-V8, vốn cho đến nay chưa được sử dụng rộng rãi để khảo sát quần xã vi khuẩn trong dạ cỏ bò đực. Việc sử dụng các cặp mồi khác nhau có thể dẫn đến sự khác biệt về độ bao phủ và khả năng định danh các đơn vị phân loại. Các vùng được sử dụng trong nghiên cứu này cũng không phải là ngoại lệ. Khi xác định số lượng sinh vật được nhận dạng ở mỗi cấp độ phân loại, ba trong số năm cấp độ phân loại được kiểm tra đã cho thấy sự phong phú hơn, chỉ ra số lượng lớn hơn các đơn vị phân loại được nhận dạng khi sử dụng đoạn khuếch đại V1-V8 để định danh 16S. Kết quả này phù hợp với các nghiên cứu trước đây, vốn đã chứng minh rằng đoạn đọc dài hơn cung cấp độ phân giải phát sinh loài cao hơn. Các đoạn khuếch đại dài hơn đã xác định được số lượng ngành lớn hơn, chẳng hạn như Proteobacteria và Verrucomicrobia. Mặc dù Proteobacteria không được coi là ngành chính trong dạ cỏ bò đực, rõ ràng là việc phân loại chỉ sử dụng vùng siêu biến V1-V3 có thể đã góp phần vào việc ước tính thấp mức độ phong phú của nó trong nghiên cứu này. Sự khác biệt về chức năng, quy trình và nền tảng tồn tại giữa hai phương pháp giải trình tự, và các yếu tố giả tạo (artifacts) từ những khác biệt đó có thể được phản ánh trong bộ dữ liệu. Việc ước tính thành phần quần xã vi khuẩn dựa trên các đoạn khuếch đại gen 16S rRNA luôn tiềm ẩn một mức độ sai lệch nhất định do mồi (primer bias). Ngoài ra, thiết bị PacBio có thể ưu tiên tải các đoạn DNA nhỏ hơn. Chúng tôi không tìm thấy bằng chứng nào cho thấy sai lệch cụ thể này ảnh hưởng đến kết quả. Dữ liệu được cắt ngắn (in silico) cho thấy một mức độ sai lệch do mồi có thể tồn tại trong dữ liệu có nguồn gốc từ các đoạn đọc PacBio, được chỉ ra bởi sự dịch chuyển trong mức độ phong phú phân loại tương đối của Proteobacteria và Bacteroidetes so với các trình tự đoạn khuếch đại V1-V8 đầy đủ. Dữ liệu hỗ trợ việc ước tính thấp Proteobacteria bằng vùng V1-V3 có nguồn gốc từ MiSeq. Dữ liệu metagenome 16S (Hình 3) đã cung cấp thêm thông tin về khả năng sai lệch do mồi và làm rõ hơn về độ phong phú tương đối của Proteobacteria và Bacteroidetes. Mặc dù có sự khác biệt (ví dụ, ước tính quá cao Proteobacteria và ước tính thấp Bacteroidetes) so với các mô tả truyền thống về quần xã vi sinh vật dạ cỏ, dữ liệu hiện tại gợi ý rằng sự thay đổi trong hồ sơ phân loại có thể không hoàn toàn do sai lệch từ mồi, mà còn liên quan đến khả năng hạn chế của đoạn khuếch đại 16S ngắn trong việc định danh chính xác các đơn vị phân loại chủ chốt. Phát hiện này có ý nghĩa quan trọng đối với các nghiên cứu về hệ vi sinh vật dạ cỏ, đặc biệt khi các nghiên cứu dựa trên DNA và RNA gần đây đã gợi ý về vai trò ngày càng tăng của Proteobacteria trong quá trình trao đổi chất ở dạ cỏ. Dữ liệu này cũng nhấn mạnh sự cần thiết phải phân tích sâu hơn và kiểm tra độ phân giải tốt hơn của quần xã vi sinh vật, không chỉ dựa vào kết quả của một đoạn 16S hoặc vùng siêu biến duy nhất. Điều này nhấn mạnh rằng việc sử dụng đồng thời nhiều bộ dữ liệu hoặc dùng một phương pháp làm đối chứng cho độ chính xác có thể cải thiện đáng kể chất lượng của các phân tích quần xã vi sinh vật dựa trên gen ribosome. Thông qua so sánh, có thể thu được một bức tranh toàn diện hơn về quần xã, từ đó hỗ trợ tốt hơn cho việc tìm hiểu các mối quan hệ và tương tác tiềm ẩn. Ngoài sự biến động về mức độ phong phú phân loại tương đối, nghiên cứu cũng kiểm tra độ chính xác phân loại giữa các vùng siêu biến được khảo sát. Thành công trong phân loại là một yếu tố cơ bản trong phân tích quần xã, nơi phần lớn dữ liệu có thể bị mất, tùy thuộc vào kỹ thuật được sử dụng trong nghiên cứu. Dữ liệu mức độ phong phú phân loại tương đối trong Hình 2 chỉ xác định các đoạn đọc không được gán tên ở cấp độ phân loại cụ thể đó; tuy nhiên, nghiên cứu cũng nhằm định lượng độ chính xác tổng thể, bao gồm các đoạn đọc không được gán tên tích lũy ở tất cả các cấp độ. Hình 4 minh họa rằng tỷ lệ các đoạn đọc 16S rRNA V1-V8 do PacBio tạo ra được phân loại thành công lớn hơn đáng kể (P < 0,05) ở mỗi cấp độ phân loại so với các đoạn đọc V1-V3 khi sử dụng Cơ sở dữ liệu gen 16S rRNA Greengenes. Điều này phù hợp với các nghiên cứu hiện tại chứng minh độ chính xác phân loại được cải thiện khi sử dụng đoạn khuếch đại 16S dài hơn. Nhìn chung, sự phong phú của các trình tự không được phân loại trong cả hai bộ dữ liệu có thể liên quan đến cơ sở dữ liệu và đặc thù của môi trường. Trong bối cảnh này, nghiên cứu đã chỉ ra rằng phân tích phát sinh loài của các vi sinh vật thường trú trong dạ cỏ có thể bị hạn chế bởi sự thiếu hụt của cơ sở dữ liệu tham chiếu. Để so sánh mức độ tương đồng giữa các quần xã được thiết lập bởi các nhóm đoạn khuếch đại thay vì giữa các mẫu, chúng tôi đã sử dụng ma trận khoảng cách UniFrac với phương pháp phân cụm UPGMA có kiểm định jackknife cũng như PCoA. Cây UPGMA (Hình 5) và biểu đồ PCoA (Hình 6) minh họa rằng có sự khác biệt giữa khoảng cách UniFrac không trọng số và có trọng số ở cả chiều dài nhánh và cách phân cụm các mẫu khi sử dụng dữ liệu trình tự từ các đoạn khuếch đại khác nhau, và những khác biệt này là có ý nghĩa thống kê. Kết quả cho thấy việc lựa chọn vùng trình tự một phần hay gần toàn bộ chiều dài của các đoạn khuếch đại được phân tích đã ảnh hưởng đáng kể đến đánh giá phát sinh loài của quần xã vi sinh vật khi sử dụng UniFrac. Những dữ liệu này chứng minh rằng các phân tích quần xã khác nhau giữa các bộ dữ liệu trình tự từ các đoạn khuếch đại khác nhau, càng củng cố thêm tầm quan trọng của việc lựa chọn cẩn thận các đoạn khuếch đại 16S trong phân tích quần xã vi sinh vật.
4.1 Kết luận
Việc lựa chọn đoạn khuếch đại mục tiêu là một bước tối quan trọng trong các nghiên cứu quần xã vi sinh vật dựa trên phân tích gen 16S rRNA. Nghiên cứu này nhấn mạnh sự khác biệt trong việc mô tả quần xã vi khuẩn dạ cỏ bò phụ thuộc vào đoạn khuếch đại 16S rRNA được sử dụng. Nghiên cứu đã chứng minh sự gia tăng về độ phân giải phát sinh loài, sự khác biệt trong phân loại định danh, độ chính xác phân loại cao hơn, cũng như sự khác biệt trong thành phần quần xã, khi sử dụng đoạn đọc 16S gần toàn bộ chiều dài trong dạ cỏ bò đực so với đoạn đọc V1-V3 ngắn hơn. Mặc dù kết quả cho thấy sự ưu việt của việc sử dụng đoạn đọc toàn bộ hoặc gần toàn bộ gen cho phân tích quần xã, các yếu tố như hạn chế về ngân sách có thể cản trở việc triển khai các nghiên cứu quy mô lớn, chuyên sâu sử dụng đoạn đọc dài từ thiết bị Pacific Biosciences RSII. Việc xây dựng các cơ sở dữ liệu vi sinh vật tham chiếu, đặc thù cho từng cộng đồng, có thể là một giải pháp thực tế hơn để nâng cao độ chính xác trong xác định hệ vi sinh vật. Trong trường hợp như vậy, các đoạn đọc dài hơn được tạo ra bởi thiết bị PacBio có thể ưu việt hơn trong việc xác định hệ vi sinh vật đặc thù cho môi trường, làm tham chiếu cho các phân tích 16S thông lượng cao, đọc ngắn hơn. Điều này cũng sẽ hỗ trợ nâng cao độ chính xác của cơ sở dữ liệu và hạn chế số lượng các đơn vị phân loại không được gán tên. Một phương pháp tiếp cận tích hợp, sử dụng đồng thời cả hai kích thước đoạn khuếch đại 16S, có thể là cần thiết để đạt được sự hiểu biết toàn diện nhất về quần xã vi khuẩn. Nghiên cứu về các hệ sinh thái phức tạp và cấu trúc phát sinh loài của quần xã vi khuẩn, như hệ vi sinh vật dạ cỏ, đã được cách mạng hóa nhờ vào sự ra đời của các nền tảng giải trình tự thông lượng cao với độ sâu lớn. Việc sử dụng các nền tảng giải trình tự thế hệ mới, chẳng hạn như MiSeq, đã cho phép tạo ra lượng lớn dữ liệu cho phân tích quần xã. Tuy nhiên, việc ứng dụng các nền tảng này thường đi kèm với việc sử dụng đoạn đọc ngắn hơn, điều này có thể cản trở việc phân tích chính xác quần xã tùy thuộc vào môi trường được kiểm tra. Nghiên cứu này đã sử dụng hai nền tảng giải trình tự để xác định hiệu quả của việc mô tả quần xã vi sinh vật dạ cỏ từ đoạn khuếch đại 16S rRNA V1-V3 ngắn hơn, với nhận định rằng nhiều nghiên cứu đã và sẽ tiếp tục sử dụng các nền tảng thông lượng cao để thu thập dữ liệu hệ vi sinh vật. Mặc dù các nền tảng này mang lại lợi thế về thời gian và khối lượng dữ liệu, nghiên cứu của chúng tôi ủng hộ việc áp dụng đa dạng các phương pháp tiếp cận ban đầu để đảm bảo tính chính xác của các phân tích quần xã về sau.
LOBI Vietnam là công ty tiên phong trong lĩnh vực Đọc trình tự gen thế hệ mới NGS (Next Generation Sequencing) và Phân tích Tin sinh học. Liên hệ hotline/Zalo 092.510.8899 để biết thêm chi tiết.
