diff --git a/src/MAT.jl b/src/MAT.jl
index a11efee..89efdb3 100644
--- a/src/MAT.jl
+++ b/src/MAT.jl
@@ -26,11 +26,12 @@ module MAT
 
 using HDF5, SparseArrays
 
+include("MAT_subsys.jl")
 include("MAT_HDF5.jl")
 include("MAT_v5.jl")
 include("MAT_v4.jl")
 
-using .MAT_HDF5, .MAT_v5, .MAT_v4
+using .MAT_HDF5, .MAT_v5, .MAT_v4, .MAT_subsys
 
 export matopen, matread, matwrite, @read, @write
 
@@ -40,7 +41,7 @@ function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Boo
     # When creating new files, create as HDF5 by default
     fs = filesize(filename)
     if cr && (tr || fs == 0)
-        return MAT_HDF5.matopen(filename, rd, wr, cr, tr, ff, compress)
+        return MAT_HDF5.matopen(filename, rd, wr, cr, tr, ff, compress, Base.ENDIAN_BOM == 0x04030201)
     elseif fs == 0
         error("File \"$filename\" does not exist and create was not specified")
     end
@@ -76,7 +77,7 @@ function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Boo
         seek(rawfid, offset)
         if read!(rawfid, Vector{UInt8}(undef, 8)) == HDF5_HEADER
             close(rawfid)
-            return MAT_HDF5.matopen(filename, rd, wr, cr, tr, ff, compress)
+            return MAT_HDF5.matopen(filename, rd, wr, cr, tr, ff, compress, endian_indicator == 0x494D)
         end
     end
 
@@ -133,6 +134,7 @@ function matread(filename::AbstractString)
     try
         vars = read(file)
     finally
+        MAT_subsys.clear_subsys!()
         close(file)
     end
     vars
@@ -165,7 +167,7 @@ function matwrite(filename::AbstractString, dict::AbstractDict{S, T}; compress::
         end
 
     else
-        
+
         file = matopen(filename, "w"; compress = compress)
         try
             for (k, v) in dict
diff --git a/src/MAT_HDF5.jl b/src/MAT_HDF5.jl
index b0d36d6..05fdede 100644
--- a/src/MAT_HDF5.jl
+++ b/src/MAT_HDF5.jl
@@ -29,6 +29,7 @@
 module MAT_HDF5
 
 using HDF5, SparseArrays
+using ..MAT_subsys
 
 import Base: names, read, write, close
 import HDF5: Reference
@@ -69,8 +70,13 @@ function close(f::MatlabHDF5File)
                 unsafe_copyto!(magicptr, idptr, length(identifier))
             end
             magic[126] = 0x02
-            magic[127] = 0x49
-            magic[128] = 0x4d
+            if Base.ENDIAN_BOM == 0x04030201
+                magic[127] = 0x49
+                magic[128] = 0x4d
+            else
+                magic[127] = 0x4d
+                magic[128] = 0x49
+            end
             rawfid = open(f.plain.filename, "r+")
             write(rawfid, magic)
             close(rawfid)
@@ -80,7 +86,7 @@ function close(f::MatlabHDF5File)
     nothing
 end
 
-function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool, compress::Bool)
+function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool, compress::Bool, endian_indicator::Bool)
     local f
     if ff && !wr
         error("Cannot append to a read-only file")
@@ -109,6 +115,11 @@ function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Boo
         fid.refcounter = length(g)-1
         close(g)
     end
+    subsys_refs = "#subsystem#"
+    if haskey(fid.plain, subsys_refs)
+        subsys_data = m_read(fid.plain[subsys_refs])
+        MAT_subsys.load_subsys!(subsys_data, endian_indicator)
+    end
     fid
 end
 
@@ -118,6 +129,7 @@ const name_type_attr_matlab = "MATLAB_class"
 const empty_attr_matlab = "MATLAB_empty"
 const sparse_attr_matlab = "MATLAB_sparse"
 const int_decode_attr_matlab = "MATLAB_int_decode"
+const object_type_attr_matlab = "MATLAB_object_decode"
 
 ### Reading
 function read_complex(dtype::HDF5.Datatype, dset::HDF5.Dataset, ::Type{T}) where T
@@ -128,6 +140,21 @@ function read_complex(dtype::HDF5.Datatype, dset::HDF5.Dataset, ::Type{T}) where
     return read(dset, Complex{T})
 end
 
+function read_cell(dset::HDF5.Dataset)
+    refs = read(dset, Reference)
+    out = Array{Any}(undef, size(refs))
+    f = HDF5.file(dset)
+    for i = 1:length(refs)
+        dset = f[refs[i]]
+        try
+            out[i] = m_read(dset)
+        finally
+            close(dset)
+        end
+    end
+    return out
+end
+
 function m_read(dset::HDF5.Dataset)
     if haskey(dset, empty_attr_matlab)
         # Empty arrays encode the dimensions as the dataset
@@ -150,36 +177,46 @@ function m_read(dset::HDF5.Dataset)
     end
 
     mattype = haskey(dset, name_type_attr_matlab) ? read_attribute(dset, name_type_attr_matlab) : "cell"
+    objecttype = haskey(dset, object_type_attr_matlab) ? read_attribute(dset, object_type_attr_matlab) : nothing
 
-    if mattype == "cell"
+    if mattype == "cell" && objecttype === nothing
         # Cell arrays, represented as an array of refs
-        refs = read(dset, Reference)
-        out = Array{Any}(undef, size(refs))
-        f = HDF5.file(dset)
-        for i = 1:length(refs)
-            dset = f[refs[i]]
-            try
-                out[i] = m_read(dset)
-            finally
-                close(dset)
-            end
+        return read_cell(dset)
+    elseif objecttype !== nothing
+        if objecttype != 3
+            @warn "MATLAB Object Type $mattype is currently not supported."
+            return missing
+        end
+        if mattype == "FileWrapper__"
+            return read_cell(dset)
+        end
+        if haskey(dset, "MATLAB_fields")
+            @warn "Enumeration Instances are not supported currently."
+            return missing
         end
-        return out
     elseif !haskey(str2type_matlab,mattype)
-        @warn "MATLAB $mattype values are currently not supported"
+        @warn "MATLAB $mattype values are currently not supported."
         return missing
     end
 
     # Regular arrays of values
     # Convert to Julia type
-    T = str2type_matlab[mattype]
+    if objecttype === nothing
+        T = str2type_matlab[mattype]
+    else
+        T = UInt32 # FIXME: Default for MATLAB objects?
+    end
 
     # Check for a COMPOUND data set, and if so handle complex numbers specially
     dtype = datatype(dset)
     try
         class_id = HDF5.API.h5t_get_class(dtype.id)
         d = class_id == HDF5.API.H5T_COMPOUND ? read_complex(dtype, dset, T) : read(dset, T)
-        length(d) == 1 ? d[1] : d
+        if objecttype !== nothing
+            return MAT_subsys.load_mcos_object(d, "MCOS")
+        else
+            return length(d) == 1 ? d[1] : d
+        end
     finally
         close(dtype)
     end
@@ -194,7 +231,11 @@ end
 
 # reading a struct, struct array, or sparse matrix
 function m_read(g::HDF5.Group)
-    mattype = read_attribute(g, name_type_attr_matlab)
+    if HDF5.name(g) == "/#subsystem#"
+        mattype = "#subsystem#"
+    else
+        mattype = read_attribute(g, name_type_attr_matlab)
+    end
     if mattype != "struct"
         # Check if this is a sparse matrix.
         fn = keys(g)
@@ -226,10 +267,11 @@ function m_read(g::HDF5.Group)
             end
             return SparseMatrixCSC(convert(Int, read_attribute(g, sparse_attr_matlab)), length(jc)-1, jc, ir, data)
         elseif mattype == "function_handle"
-            @warn "MATLAB $mattype values are currently not supported"
-            return missing
+            # Fall through
         else
-            @warn "Unknown non-struct group of type $mattype detected; attempting to read as struct"
+            if mattype != "#subsystem#"
+                @warn "Unknown non-struct group of type $mattype detected; attempting to read as struct"
+            end
         end
     end
     if haskey(g, "MATLAB_fields")
diff --git a/src/MAT_subsys.jl b/src/MAT_subsys.jl
new file mode 100644
index 0000000..7d7eef2
--- /dev/null
+++ b/src/MAT_subsys.jl
@@ -0,0 +1,297 @@
+# MAT_subsys.jl
+# Tools for processing MAT-file subsystem data in Julia
+#
+# Copyright (C) 2025    Nithin Lakshmisha
+#
+# Permission is hereby granted, free of charge, to any person obtaining
+# a copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+#
+# The above copyright notice and this permission notice shall be
+# included in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+module MAT_subsys
+
+const FWRAP_VERSION = 4
+
+mutable struct Subsys
+    num_names::UInt32
+    mcos_names::Vector{String}
+    class_id_metadata::Vector{UInt32}
+    object_id_metadata::Vector{UInt32}
+    saveobj_prop_metadata::Vector{UInt32}
+    obj_prop_metadata::Vector{UInt32}
+    dynprop_metadata::Vector{UInt32}
+    _u6_metadata::Vector{UInt32}
+    _u7_metadata::Vector{UInt32}
+    prop_vals_saved::Vector{Any}
+    _c3::Any
+    _c2::Any
+    prop_vals_defaults::Any
+    handle_data::Any
+    java_data::Any
+
+    Subsys() = new(
+        UInt32(0),
+        String[],
+        UInt32[],
+        UInt32[],
+        UInt32[],
+        UInt32[],
+        UInt32[],
+        UInt32[],
+        UInt32[],
+        Any[],
+        nothing,
+        nothing,
+        nothing,
+        nothing,
+        nothing
+    )
+end
+
+const subsys_cache = Ref{Union{Nothing,Subsys}}(nothing)
+const object_cache = Ref{Union{Nothing, Dict{UInt32, Dict{String,Any}}}}(nothing)
+
+function clear_subsys!()
+    subsys_cache[] = nothing
+    object_cache[] = nothing
+end
+
+function load_subsys!(subsystem_data::Dict{String,Any}, swap_bytes::Bool)
+    subsys_cache[] = Subsys()
+    object_cache[] = Dict{UInt32, Dict{String,Any}}()
+    subsys_cache[].handle_data = get(subsystem_data, "handle", nothing)
+    subsys_cache[].java_data = get(subsystem_data, "java", nothing)
+    mcos_data = get(subsystem_data, "MCOS", nothing)
+    if mcos_data === nothing
+        return
+    end
+
+    if mcos_data isa Tuple
+        # Backward compatibility with MAT_v5
+        mcos_data = mcos_data[2]
+    end
+    fwrap_metadata = vec(mcos_data[1, 1])
+
+    # FIXME: Is this the best way to read?
+    # Integers are written as uint8 (with swap), interpret as uint32
+    version = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[1:4]) : fwrap_metadata[1:4])[1]
+    if version <= 1 || version > FWRAP_VERSION
+        error("Cannot read subsystem: Unsupported FileWrapper version: $version")
+    end
+
+    subsys_cache[].num_names = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[5:8]) : fwrap_metadata[5:8])[1]
+    region_offsets = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[9:40]) : fwrap_metadata[9:40])
+
+    # Class and Property Names stored as list of null-terminated strings
+    start = 41
+    pos = start
+    name_count = 0
+    while name_count < subsys_cache[].num_names
+        if fwrap_metadata[pos] == 0x00
+            push!(subsys_cache[].mcos_names, String(fwrap_metadata[start:pos-1]))
+            name_count += 1
+            start = pos + 1
+            if name_count == subsys_cache[].num_names
+                break
+            end
+        end
+        pos += 1
+    end
+
+    subsys_cache[].class_id_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[1]+1:region_offsets[2]]) : fwrap_metadata[region_offsets[1]+1:region_offsets[2]])
+    subsys_cache[].saveobj_prop_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[2]+1:region_offsets[3]]) : fwrap_metadata[region_offsets[2]+1:region_offsets[3]])
+    subsys_cache[].object_id_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[3]+1:region_offsets[4]]) : fwrap_metadata[region_offsets[3]+1:region_offsets[4]])
+    subsys_cache[].obj_prop_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[4]+1:region_offsets[5]]) : fwrap_metadata[region_offsets[4]+1:region_offsets[5]])
+    subsys_cache[].dynprop_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[5]+1:region_offsets[6]]) : fwrap_metadata[region_offsets[5]+1:region_offsets[6]])
+
+    if region_offsets[7] != 0
+        subsys_cache[]._u6_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[6]+1:region_offsets[7]]) : fwrap_metadata[region_offsets[6]+1:region_offsets[7]])
+    end
+
+    if region_offsets[8] != 0
+        subsys_cache[]._u7_metadata = reinterpret(UInt32, swap_bytes ? reverse(fwrap_metadata[region_offsets[7]+1:region_offsets[8]]) : fwrap_metadata[region_offsets[7]+1:region_offsets[8]])
+    end
+
+    if version == 2
+        subsys_cache[].prop_vals_saved = mcos_data[3:end-1, 1]
+    elseif version == 3
+        subsys_cache[].prop_vals_saved = mcos_data[3:end-2, 1]
+        subsys_cache[]._c2 = mcos_data[end-1, 1]
+    else
+        subsys_cache[].prop_vals_saved = mcos_data[3:end-3, 1]
+        subsys_cache[]._c3 = mcos_data[end-2, 1]
+    end
+
+    subsys_cache[].prop_vals_defaults = mcos_data[end, 1]
+end
+
+function get_classname(class_id::UInt32)
+    namespace_idx = subsys_cache[].class_id_metadata[class_id*4+1]
+    classname_idx = subsys_cache[].class_id_metadata[class_id*4+2]
+
+    namespace = if namespace_idx == 0
+        ""
+    else
+        subsys_cache[].mcos_names[namespace_idx] * "."
+    end
+
+    classname = namespace * subsys_cache[].mcos_names[classname_idx]
+    return classname
+end
+
+function get_object_metadata(object_id::UInt32)
+    return subsys_cache[].object_id_metadata[object_id*6+1:object_id*6+6]
+end
+
+function get_default_properties(class_id::UInt32)
+    # FIXME Should we use deepcopy here
+    return copy(subsys_cache[].prop_vals_defaults[class_id+1, 1])
+end
+
+function get_property_idxs(obj_type_id::UInt32, saveobj_ret_type::Bool)
+    prop_field_idxs = saveobj_ret_type ? subsys_cache[].saveobj_prop_metadata : subsys_cache[].obj_prop_metadata
+    nfields = 3
+    offset = 1
+    while obj_type_id > 0
+        nprops = prop_field_idxs[offset]
+        offset += 1 + (nfields * nprops)
+        offset += (offset + 1) % 2  # Padding
+        obj_type_id -= 1
+    end
+    nprops = prop_field_idxs[offset]
+    offset += 1
+    return prop_field_idxs[offset:offset+nprops*nfields-1]
+end
+
+function find_nested_prop(prop_value::Any)
+    # Hacky way to find a nested object
+    # Nested objects are stored as a uint32 Matrix with a unique signature
+    # MATLAB probably uses some kind of placeholders to decode
+    # But this should work here
+    if prop_value isa Dict
+        # Handle nested objects in a dictionary (struct)
+        for (key, value) in prop_value
+            prop_value[key] = find_nested_prop(value)
+        end
+    end
+
+    if prop_value isa Matrix{Any}
+        # Handle nested objects in a Cell
+        for i in eachindex(prop_value)
+            prop_value[i] = find_nested_prop(prop_value[i])
+        end
+    end
+
+    if prop_value isa Matrix{UInt32} && prop_value[1,1] == 0xdd000000
+        # MATLAB identifies any uint32 array with first value 0xdd000000 as an MCOS object
+        return load_mcos_object(prop_value, "MCOS")
+    end
+
+    return prop_value
+end
+
+function get_saved_properties(obj_type_id::UInt32, saveobj_ret_type::Bool)
+    save_prop_map = Dict{String,Any}()
+    prop_field_idxs = get_property_idxs(obj_type_id, saveobj_ret_type)
+    nprops = length(prop_field_idxs) ÷ 3
+    for i in 0:nprops-1
+        prop_name = subsys_cache[].mcos_names[prop_field_idxs[i*3+1]]
+        prop_type = prop_field_idxs[i*3+2]
+        if prop_type == 0
+            prop_value = subsys_cache[].mcos_names[prop_field_idxs[i*3+3]]
+        elseif prop_type == 1
+            prop_value = subsys_cache[].prop_vals_saved[prop_field_idxs[i*3+3]+1]
+        elseif prop_type == 2
+            prop_value = prop_field_idxs[i*3+3]
+        else
+            error("Unknown property type ID: $prop_type encountered during deserialization")
+        end
+        save_prop_map[prop_name] = find_nested_prop(prop_value)
+    end
+    return save_prop_map
+end
+
+function get_properties(object_id::UInt32)
+    if object_id == 0
+        return Dict{String,Any}()
+    end
+
+    class_id, _, _, saveobj_id, normobj_id, _ = get_object_metadata(object_id)
+    if saveobj_id != 0
+        saveobj_ret_type = true
+        obj_type_id = saveobj_id
+    else
+        saveobj_ret_type = false
+        obj_type_id = normobj_id
+    end
+
+    prop_map = get_default_properties(class_id)
+    merge!(prop_map, get_saved_properties(obj_type_id, saveobj_ret_type))
+    # TODO: Add dynamic properties
+    return prop_map
+end
+
+function load_mcos_object(metadata::Any, type_name::String)
+    if type_name != "MCOS"
+        @warn "Loading Type:$type_name is not implemented. Returning metadata."
+        return metadata
+    end
+
+    if isa(metadata, Dict)
+        # TODO: Load Enumeration Instances
+        @warn "Loading enumeration instances are not supported. Returning Metadata"
+        return metadata
+    end
+
+    if !(metadata isa Array{UInt32})
+        @warn "Expected MCOS metadata to be an Array{UInt32}, got $(typeof(metadata)). Returning metadata."
+        return metadata
+    end
+
+    if metadata[1, 1] != 0xDD000000
+        @warn "MCOS object metadata is corrupted. Returning raw data."
+        return metadata
+    end
+
+    ndims = metadata[2, 1]
+    dims = metadata[3:2+ndims, 1]
+    nobjects = prod(dims)
+    object_ids = metadata[3+ndims:2+ndims+nobjects, 1]
+
+    class_id = metadata[end, 1]
+    classname = get_classname(class_id)
+
+    object_arr = Array{Dict{String,Any}}(undef, convert(Vector{Int}, dims)...)
+
+    for i = 1:length(object_arr)
+        oid = object_ids[i]
+        if haskey(object_cache[], oid)
+            prop_dict = object_cache[][oid]
+        else
+            prop_dict = Dict{String,Any}()
+            object_cache[][oid] = prop_dict
+            merge!(prop_dict, get_properties(oid))
+            prop_dict["__class__"] = classname
+        end
+        object_arr[i] = prop_dict
+    end
+
+    return object_arr
+
+end
+
+end
\ No newline at end of file
diff --git a/src/MAT_v5.jl b/src/MAT_v5.jl
index 7370dd4..c412a71 100644
--- a/src/MAT_v5.jl
+++ b/src/MAT_v5.jl
@@ -29,6 +29,8 @@ module MAT_v5
 using CodecZlib, BufferedStreams, HDF5, SparseArrays
 import Base: read, write, close
 
+using ..MAT_subsys
+
 round_uint8(data) = round.(UInt8, data)
 complex_array(a, b) = complex.(a, b)
 
@@ -246,7 +248,7 @@ function read_sparse(f::IO, swap_bytes::Bool, dimensions::Vector{Int32}, flags::
     end
     if length(ir) > length(pr)
         # Fix for Issue #169, xref https://github.com/JuliaLang/julia/pull/40523
-        #= 
+        #=
         # The following expression must be obeyed according to
         # https://github.com/JuliaLang/julia/blob/b3e4341d43da32f4ab6087230d98d00b89c8c004/stdlib/SparseArrays/src/sparsematrix.jl#L86-L90
         @debug "SparseMatrixCSC" m n jc ir pr
@@ -311,6 +313,18 @@ function read_string(f::IO, swap_bytes::Bool, dimensions::Vector{Int32})
     data
 end
 
+function read_opaque(f::IO, swap_bytes::Bool)
+    type_name = String(read_element(f, swap_bytes, UInt8))
+    classname = String(read_element(f, swap_bytes, UInt8))
+
+    if classname == "FileWrapper__"
+        return read_matrix(f, swap_bytes)
+    end
+
+    _, metadata = read_matrix(f, swap_bytes)
+    return MAT_subsys.load_mcos_object(metadata, type_name)
+end
+
 # Read matrix data
 function read_matrix(f::IO, swap_bytes::Bool)
     (dtype, nbytes) = read_header(f, swap_bytes)
@@ -332,15 +346,10 @@ function read_matrix(f::IO, swap_bytes::Bool)
     flags = read_element(f, swap_bytes, UInt32)
     class = flags[1] & 0xFF
 
-    if class == mxOPAQUE_CLASS
-        s0 = read_data(f, swap_bytes)
-        s1 = read_data(f, swap_bytes)
-        s2 = read_data(f, swap_bytes)
-        arr = read_matrix(f, swap_bytes)
-        return ("__opaque__", Dict("s0"=>s0, "s1"=>s1, "s2"=>s2, "arr"=>arr))
+    if class != mxOPAQUE_CLASS
+        dimensions = read_element(f, swap_bytes, Int32)
     end
 
-    dimensions = read_element(f, swap_bytes, Int32)
     name = String(read_element(f, swap_bytes, UInt8))
 
     local data
@@ -354,6 +363,8 @@ function read_matrix(f::IO, swap_bytes::Bool)
         data = read_string(f, swap_bytes, dimensions)
     elseif class == mxFUNCTION_CLASS
         data = read_matrix(f, swap_bytes)
+    elseif class == mxOPAQUE_CLASS
+        data = read_opaque(f, swap_bytes)
     else
         if (flags[1] & (1 << 9)) != 0 # logical
             data = read_data(f, swap_bytes, Bool, dimensions)
@@ -375,14 +386,41 @@ matopen(ios::IOStream, endian_indicator::UInt16) =
 
 # Read whole MAT file
 function read(matfile::Matlabv5File)
-    seek(matfile.ios, 128)
     vars = Dict{String, Any}()
+
+    seek(matfile.ios, 116)
+    subsys_offset = read_bswap(matfile.ios, matfile.swap_bytes, UInt64)
+    if subsys_offset == 0x2020202020202020
+        subsys_offset = 0
+    end
+    if subsys_offset != 0
+        read_subsystem(matfile.ios, matfile.swap_bytes, subsys_offset)
+    end
+
+    seek(matfile.ios, 128)
     while !eof(matfile.ios)
+        offset = position(matfile.ios)
+        if offset == subsys_offset
+            # Skip reading subsystem again
+            (_, nbytes) = read_header(matfile.ios, matfile.swap_bytes)
+            skip(matfile.ios, nbytes)
+            continue
+        end
         (name, data) = read_matrix(matfile.ios, matfile.swap_bytes)
         vars[name] = data
     end
     vars
 end
+
+function read_subsystem(ios::IOStream, swap_bytes::Bool, offset::UInt64)
+    seek(ios, offset)
+    (_, subsystem_data) = read_matrix(ios, swap_bytes)
+    buf = IOBuffer(vec(subsystem_data))
+    seek(buf, 8) # Skip subsystem header
+    _, subsys_data = read_matrix(buf, swap_bytes)
+    MAT_subsys.load_subsys!(subsys_data, swap_bytes)
+end
+
 # Read only variable names from an HDF5 file
 function getvarnames(matfile::Matlabv5File)
     if !isdefined(matfile, :varnames)
diff --git a/test/read.jl b/test/read.jl
index 4d8c9d4..64b2761 100644
--- a/test/read.jl
+++ b/test/read.jl
@@ -214,21 +214,46 @@ let objtestfile = "figure.fig"
 end
 
 # test reading file containing Matlab function handle, table, and datetime objects
-# since we don't support these objects, just make sure that there are no errors
-# reading the file and that the variables are there and replaced with `missing`
 let objtestfile = "function_handles.mat"
     vars = matread(joinpath(dirname(@__FILE__), "v7.3", objtestfile))
     @test "sin" in keys(vars)
-    @test ismissing(vars["sin"])
+    @test typeof(vars["sin"]) == Dict{String, Any}
+    @test Set(keys(vars["sin"])) == Set(["function_handle", "sentinel", "separator", "matlabroot"])
     @test "anonymous" in keys(vars)
-    @test ismissing(vars["anonymous"])
+    @test typeof(vars["anonymous"]) == Dict{String, Any}
+    @test Set(keys(vars["anonymous"])) == Set(["function_handle", "sentinel", "separator", "matlabroot"])
 end
-let objtestfile = "struct_table_datetime.mat"
-    vars = matread(joinpath(dirname(@__FILE__), "v7.3", objtestfile))["s"]
-    @test "testTable" in keys(vars)
-    @test ismissing(vars["testTable"])
-    @test "testDatetime" in keys(vars)
-    @test ismissing(vars["testDatetime"])
+
+for format in ["v7", "v7.3"]
+    let objtestfile = "struct_table_datetime.mat"
+        vars = matread(joinpath(dirname(@__FILE__), format, objtestfile))["s"]
+        @test "testTable" in keys(vars)
+        @test size(vars["testTable"]) == (1, 1)
+        @test Set(keys(vars["testTable"][1, 1])) == Set(["__class__", "props", "varnames", "nrows", "data", "rownames", "ndims", "nvars"])
+        @test vars["testTable"][1, 1]["__class__"] == "table"
+        @test vars["testTable"][1, 1]["ndims"] === 2.0
+        @test vars["testTable"][1, 1]["nvars"] === 5.0
+        @test vars["testTable"][1, 1]["nrows"] === 3.0
+        @test vars["testTable"][1, 1]["data"][1, 1] == reshape([1261.0, 547.0, 3489.0], 3, 1)
+        @test vars["testTable"][1, 1]["data"][1, 2][1, 1]["__class__"] == "string"
+        @test vars["testTable"][1, 1]["data"][1, 3][1, 1]["__class__"] == "datetime"
+        @test vars["testTable"][1, 1]["data"][1, 4][1, 1]["__class__"] == "categorical"
+        @test vars["testTable"][1, 1]["data"][1, 5][1, 1]["__class__"] == "string"
+
+        @test "testDatetime" in keys(vars)
+        @test size(vars["testDatetime"]) == (1, 1)
+        if format == "v7.3"
+            @test Set(keys(vars["testDatetime"][1, 1])) == Set(["__class__", "tz", "data", "fmt", "isDateOnly"])
+            @test vars["testDatetime"][1, 1]["isDateOnly"] === false
+        else
+            # MATLAB removed property "isDateOnly" in later versions
+            @test Set(keys(vars["testDatetime"][1, 1])) == Set(["__class__", "tz", "data", "fmt"])
+        end
+        @test vars["testDatetime"][1, 1]["__class__"] == "datetime"
+        @test vars["testDatetime"][1, 1]["tz"] === ""
+        @test vars["testDatetime"][1, 1]["fmt"] === ""
+        @test vars["testDatetime"][1, 1]["data"] === 1.575304969634e12 + 0.0im
+    end
 end
 
 # test reading of old-style Matlab object in v7.3 format
diff --git a/test/v7/struct_table_datetime.mat b/test/v7/struct_table_datetime.mat
new file mode 100644
index 0000000..6a4885a
Binary files /dev/null and b/test/v7/struct_table_datetime.mat differ