How to test differences (over time and between treatments) of a specific species in DNA metabarcoding sequencing data? I have DNA metabarcoding sequencing data in the following format:




plot
Time_point
reads_species_A
Reads_species_B
reads_species_C




1
T1
0
245
65


2
T1
48
455
0


3
T1
15
5
10


1
T3
153
23
564


2
T3
448
468
48


3
T3
753
47
5




The read counts are DNA sequencing read counts and the dataset is rarefied to have an equal number of reads per sample. Sequencing data is compositional in nature. Following article describes this compositional nature as follows: "There is increasing awareness that microbiome datasets generated by HTS are compositional because they have an arbitrary total imposed by the instrument. " https://www.frontiersin.org/articles/10.3389/fmicb.2017.02224/full
For the actual dataset, we have 66 plots which are sampled at time point "T1" and time point "T3". They are thus paired. There are also more treatments in the actual dataset, which I omitted for simplicity.
For now, I'm specifically interested in species A, and I want to test if the read count of this species differs between the two time points. As you could see it as count data (counts of reads), I know a poisson regression could be used for this. However, is this the right way to go here?
For count/frequency data, also Chi-square tests are used, but I'm not sure if it is appropriate here. I'm not really sure when to use it. Would a Chi-Square Goodness of Fit test be superior in this context?
Or would a paired Wilcoxon signed rank test be appropriate for this problem? Or logtransform the data, and see if the assumptions for a e.g. a paired t test are met?
I know that differential abundance tests (e.g. ancom, ALDeX2,..) can be used to see if there are taxa more abundant in a certain treatment than in another treatment.
However, now I'm thus looking for a way to test if a specific species, is different between time points/treatments. I also say treatments here, because for another question, I'd like to test differences of a single species between treatments (treatments no shown).
One could do a differential abundance test, and see if there is a difference in the particular species, but I'd rather select the data and they perform an appropriate test.
 A: There typically is greater dispersion in sequencing count data than would be expected from a Poisson distribution, where the variance must equal the mean. Negative-binomial regression often is used to evaluate such data, as in the Bioconductor edgeR and DESeq2 packages.
It will probably be most efficient to adapt the tools in a package like one of those to analyze all of your data together.* They can accept complex experimental designs, although you might have to do some work to get the models specified properly. That uses all of your data together, which is typically more powerful than a set of separate tests on species by species. If there's one particular species of primary interest, however, you certainly can examine that with a standard negative-binomial regression.
You will have to account for the multiple measurements within plots. In principle you could treat the plots as separate fixed effects, but in your case that leads to a lot of extra coefficients to estimate. The limma package allows you to specify things like plots as random effects for count data, greatly decreasing the number of coefficients.
With respect to compositional-data in metabarcoding, that might not be something you need to worry about in your application. If all sequencing reads could be mapped to a species, then you could only independently evaluate 1 fewer than the total number of species. In that case, the results on the last species would be completely determined by the results on all the others. I suspect that many reads aren't successfully mapped, however, so that you aren't stuck with that limitation.

*Try looking for packages specific to analyzing metabarcoding data on Bioconductor.
