Multiple, Binomial Dependent Variables for GLM (or LME4) in R

Question

everyone. I'm pretty new to R. I've been trying to educate myself about this issue, but I've continued to run into road blocks.

I have a data set with two categorical, independent variables (habitat (1,2,3) and site (1,2,3,4,5). My response variables are the presence or absence of AFLP loci. I have 96 loci, and I want to determine which, if any, of these loci are significantly associated with habitat (site is a random effect). Each of the loci can be assumed to be independent from the others.

As far as relevancy to other researchers, this should be a problem that people trying to analyze molecular data with GLM or LME will begin to run into more.

Here is my code:

##Independent variables
Site=AFLP$Site ##AFLP is my data file
Habitat=AFLP$Habitat
##Dependent variable
Loci=AFLP[,4:99]

##Establishing matrix of variables
mydata <- cbind(Site, Habitat, Loci)

##glm
model1 <- glm(Loci ~ (1|Site)+Habitat, data=mydata, family="binomial")

I get this error: Error in model.frame.default(formula = Loci ~ (1 | Site) + Habitat, data = mydata, : invalid type (list) for variable 'Loci'

I know this error is associated with the data type of Loci; however, I've tried a bunch of things and still can't figure out how to correctly address the issue.

My problem seems to be similar to the ones in the below links, but again, I haven't been able to figure out how to apply this information to my data set.

http://stackoverflow.com/questions/18067519/using-r-to-do-a-regression-with-multiple-dependent-and-multiple-independent-vari

https://stats.stackexchange.com/questions/26585/how-to-do-a-generalized-linear-model-with-multiple-dependent-variables-in-r

Thank you in advance. If this turns out to have a simple answer, I apologize for taking up space. I have been Googling and trying to educate myself, and I haven't made any head-way.