You need to Iterate in ModelBuilder.
This will allow you to cycle through each species and process that before moving on to the next.
There are several Iterate options, such as Iterate Datasets or Iterate Feature Classes.
I used a list of values (polygon of species) to perform the iterate and I solved the problem of the overwriting output insering the iteration value (%i%) in the name of the output generated by spatial join (i.e output%i%). Then I used the Join Field tool to construct the final matrix.
I am also very interested in your approach. However, you said that you "created a model (see the attached file) with "a list of values" in the species' parameter to perform the model to each species."
However, I can't access the attachment. Could you please make the model available so that i can reproduce it?
I have a set of polygons and I am interested in creating a presence/absence matrix for the polygons using grid cells generated by fishnet.
Thanks and kind regards