Dear all, 

In our study, we examined fungal communities using high-throughput amplicon sequencing (HTAS) of internal transcribed spacer 2 (ITS2) region in extracted total RNA from environmental samples (plants). We had three different treatments and performed an analysis of variance (ANOVA) to investigate which OTU's significantly different in abundance among experimental factors after Bonferroni correction.

We received the comments from the reviewers:

R1. ANOVA is not the most appropriate test for identifying OTUS that differ significantly between experimental factors, as there are issues with the effects of uneven sequencing depth and subsequent rarefactions  Instead consider implementing some of the differential expression analyses in the DESeq2 package for R.

R2. The DeSeq2 method should be used to highlight differentially abundant OTUs.

Which method would you use to test the effects of experimental factors on fungal community composition? Do you know any published example of this? 

Thanks a lot

DESeq2 can be used to compare counts across groups of samples or other complex experimental designs. Take a look at the paper and the software vignette. The input is a count matrix, and an Negative Binomial model is fit per row (OTU in your case), which takes into account a normalizing offset for sequencing depth. There is another Bioconductor package, phyloseq, which I believe has some DESeq2 workflows for OTU count tables that you can take a look at.