Using Multiple Bounding Latitudes, Not Boxes, to Clip a Polygon Shapefile in R

Question

I am working on a project in R that requires me to divide about 8 shapefile polygons with yearly data into 15 regions, each separated by latitude. Since I'm working only in shallow waters, it doesn't matter to me the cutoffs for the longitude, just latitude. I then want to look at how the area changes year to year.

I am trying to figure out code to segment each of these years' shapefiles by using 15 bounding boxes.

Problem 1: I have found code for creating a new shapefile from one bounding box, but not by several.

Problem 2: I am not sure bounding boxes are the right way to go- is it possible to use bounding latitudes only?

In order to do this I first attempted to make bounding boxes but I have a problem: even though the both the polygon projection and map projection is correct, the boxes appear as straight lines on the map, not following the latitudinal curve of the map (this is possibly because I chose only one point as the latitude, not multiple samples?). See attached photo of boxes imaged in ArcGIS.Though this would provide a rough estimate, would it end up cutting out certain areas, or adding in certain area that ultimately I won't want?

I am wondering if there is an alternative way to segment these shapefiles into new polygons using R and latitudinal cutoffs only instead of bounding boxes? Essentially there will be 15 breakpoints that would segment my polygon (ex: polygon one will be only the shapefile area found between 933707.9m N and 1038634.6m N) into regions that don't require a longitude argument. If this is the case, would it solve my problem with the latitude curve? Or, does it not matter that it appears to look this way on Arc, but will still function as a bounding box?

Here's how I created the bounding boxes, which I did pulling lat/long from a .CSV file table and in a loop:

library(sp)
library(rgdal)
library(raster)
for(i in 1:nrow(KelpBounds)){
  
  XCoords <- c(KelpBounds[i,6],KelpBounds[i,8], KelpBounds[i,9], KelpBounds [i,7])
  XCoords
  
  YCoords <- c(KelpBounds[i,4],KelpBounds[i,5], KelpBounds[i,2], KelpBounds [i,3])
  YCoords
  
  xym <- cbind(XCoords, YCoords)
  xym
  
  
  p = Polygon(xym)
  ps = Polygons(list(p),1)
  sps = SpatialPolygons(list(ps))
  plot(sps)
  
  proj4string(sps) = CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0
                       +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
  
  proj4string(sps)
  
  data = data.frame(f=99.9)
  spdf = SpatialPolygonsDataFrame(sps,data)
  spdf
  summary(spdf)
  
 ## writeOGR(spdf, dsn='NewRegions', layer= paste0("Region_",i,"_Matrix") 
       ##   driver= "ESRI Shapefile") 
  
  dsn <- layer <- gsub(".csv","",i)
  writeOGR(spdf, dsn, layer, driver="ESRI Shapefile")
  
}

thank you!