Yes, that's exactly right! Take a look at the following example -- first we take your species data frame and convert it to an sf object by using st_as_sf() and specifying the column names of the coordinates and the coordinate reference system. Then, st_distance() creates a matrix of the distances between all the points.
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
species <- structure(list(longitude = c(-106.425668, -96.111617, -96.5,
-96.5, -96.5, -96.5), latitude = c(23.18127, 15.779873, 19.68,
19.68, 19.68, 19.68)), row.names = c(NA, 6L), class = "data.frame")
species_sf <- st_as_sf(species, coords = c("longitude", "latitude"), crs = 4326)
st_distance(species_sf)
#> Units: [m]
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 1356919.2 1099211.9 1099211.9 1099211.9 1099211.9
#> [2,] 1356919 0.0 433617.2 433617.2 433617.2 433617.2
#> [3,] 1099212 433617.2 0.0 0.0 0.0 0.0
#> [4,] 1099212 433617.2 0.0 0.0 0.0 0.0
#> [5,] 1099212 433617.2 0.0 0.0 0.0 0.0
#> [6,] 1099212 433617.2 0.0 0.0 0.0 0.0
Created on 2020-06-11 by the reprex package (v0.3.0)